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SECTION  I 
INTRODUCTION 


The  FORTRAN  IV  program  0/^WNA  solves  the  set  of  partial  differential  equations 
governing  the  flow  on  the  axis  of  symmetry  between  the  blast  wave  and  the  Mach 
disc  of  a muzzle  blast  field.  A complete  description  of  the  method  of  solution 
of  the  governing  equations  and  statement  of  the  boundary  conditions  is  given 
in  Reference  ( 1 ) . Program  DAWNA  is  a revised  edition  of  program  DAWN  pre- 
viously described  in  Reference  (2).  The  two  main  refinements  included  in 
program  DAWNA  are  the  following: 

1)  An  initialization  technique  has  been  developed  that  allows  the 
computer  program  to  self-start  without  the  use  of  empirical 
relationships  (such  as  the  initial  locations  of  the  flow  dis- 
continuities). 

2)  An  acoustic  analysis  which  enables  the  continuation  of  the  blast 
field  calculation  to  very  late  times. 

Besides  the  major  analytical  refinements  presented  above,  various  other  mod- 
ifications have  been  included  in  the  present  computer  code  which  minimize 
possible  souices  of  numerical  error  as  well  as  making  the  program  more  con- 
venient for  ihe  user.  The  more  pertinent  of  these  include  the  use  of  auto- 
matic grid  subdivision,  a streamline  trace  for  the  determination  of  the  pro- 
perties upstream  of  the  Mach  disc  and  the  option  of  choosing  the  units  of 
both  the  input  and  output  parameters. 
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SECTION  II 
PROGRAM  STRUCTURE 


The  main  program  acts  an  as  executive  routine  which  calls  the  principal  sub- 
routines, as  follows.  Subroutine  INDATA  Is  called  once  to  define  the  finite 
difference  grid  and  initialize  the  dependent  variables.  Subroutines  SHOCK, 
CONTACT  and  INT  PT  are  called  sequentially  to  advance  the  solution  from  time 
t to  time  t+At.  The  dependent  variables  are  then  reinitialized,  printed  and 
the  above  sequence  repeated.  This  loop  is  terminated  when  the  selected  number 
of  time  steps  have  been  completed,  at  which  time  the  final  solution  can  be 
written  out  in  TAPE12  (at  the  option  of  the  user)  and  the  program  execution 
is  stopped. 

Subroutines  SHOCK  and  CONTACT  accomplish  the  solution  at  the  two  shocks  and  the 
contact  surface  by  the  method  of  characteristic  technique  described  in  Refer- 
ence (1).  The  three  surfaces  of  discontinuity  are  identified  by  the  indices 
1,  2 and  3 which  refer  to  the  Mach  disc,  contact  surface  and  blast  wave,  re- 
spectively. Subroutine  INT  PT  accomplishes  the  solution  at  the  interior  grid 
points  by  the  finite  difference  algorithm  devised  by  MacCormack  (Reference  3). 
Within  subroutine  INT  PT  ths  index  K«1  denotes  the  solution  between  the  Mach 
disc  and  contact  and  K-2  denotes  the  solution  between  the  contact  and  the 
blast  wave.  The  indices  L00P-O  and  1 refer  to  the  first  and  second  iterates 
of  the  MacCormack  algorithm. 

The  subroutines  associated  with  the  acoustic  wave  computation  are  called  from 
subroutine  AC0US  which  only  returns  control  to  the  main  program  at  the  com- 
pletion of  the  run.  • 

The  functions  of  all  the  subroutines  are  summarized  in  Table  I.  The  main 
program  variables  are  identified  in  Table  II. 
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NAME 

1. 

2. 

3. 

k. 

5. 

6. 

7. 
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TABLE  I 

SUBROUTINES  AND  FUNCTIONS 
DESCRIPTION 


CL 

Calculates  flow  properties  on  the  centerline  of  the  super- 
sonic Jet  exhaust  plume  from  the  muzzle  exit  conditions 
and  the  centerline  Hach  number  at  the  location  of  the  Mach 
disc. 

CONTACT 

Calculates  the  position  and  velocity  of  the  contact  sur- 
face and  the  flow  properties  on  both  sides  of  the  contact 
surface. 

FS 

Sets  the  flow  variables  on  the  upstream  side  of  the  blast 
wave  to  the  specified  ambient  (I.e.,  "infinity")  condi- 
tions. 

INDATA 

Reads  the  input  data  and  defines  accordingly  the  finite  dlf 
ference  grid  and  the  initial  values  of  the  dependent  vari- 
ables. 

INT  PT 

Computes  the  new  finite  difference  solution  at  the  interior 
grid  points  and  the  new  location  of  the  grid  points  at  each 
succeeding  time  step. 

LP0INT 

Locates  the  intersection  of  a characteristic  surface  and  a 
time  plane,  and  performs  the  necessary  interpolations  of 
variables  from  the  solution  at  the  grid  points. 

NUZZLE 

Determines  the  pressure,  speed  of  sound  and  Hach  number 
at  the  muzzle  exit  as  a function  of  time  from  a spline-fit 
of  the  corresponding  tabular  input  data.  ■ 
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NAME 

8. 

9. 

10. 

11. 

12. 

13. 

1*. 

15. 

16. 
17. 


DESCRIPTION 


PUNCH  Writes  the  final  time  step  on  TAPE12  for  restarting  the 

calculation  (see  Section  1 1 1(B)  for  an  explanation  of 
the  restart  capability  of  the  program). 

SETN  Reinitializes  the  arrays  in  which  the  new  solution  at 

time  t+At  will  be  stored. 

SHlCK  Calculates  the  position  and  velocity  of  a moving  shock 

and  the  flow  properties  on  the  downstream  side  of  the 
shock.  The  index  K in  the  calling  sequence  is  used  to 
denote  the  Mach  disc  (K-l)  and  the  blast  wave  (K<-3). 

SPLINE  Performs  a spline-fit  of  tabular  data  (see  Reference  *) . 

SPLINT  Uses  the  spline-fit  coefficients  to  interpolate  data  at 

arbitrary  values  of  the  independent  argument.  (First 
and  second  derivatives  of  the  dependent  variable  are  also 
calculated,  but  not  used  in  the  present  program.) 

STEP  Computes  the  maximum  allowable  time  step,  At,  based  on 

the  Courant,  Friedrichs  and  Lewy  criterion. 

INTER  Interpolates  data  when  adding  or  deleting  a grid  point. 

INIT  Determines  the  initial  location  and  fluid  properties  of 

the  Mach  disc,  contact  surface  and  blast  wave. 

AC0US  Calls  the  various  subroutines  associated  with  the  integra- 

tion of  the  acoustic  wave  equation. 

INTI A Performs  the  integration  of  the  acoustic  wave  equation  using 

a fourth-order  Runke-Kutta  method. 
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NAME  DESCRIPTION 

18.  DERSUB  Subroutine  which  evaluates  the  derivative  used  In  the 

Runge-Kutta  integration. 

19.  CHSUB  Dummy  subroutine  called  from  subroutine  I NT 1 A . 


FIRTRAN  NAME 
IR(!) 


IR(2) 


Z 

P 

U 

s 

A 

RH 

W 

GAM1 , GAM2 


TABLE  II 

PRINCIPAL  PROGRAM  VARIABLES 
DEFINITION 

Number  of  grid  points  In  region  1 (l.e.,  from  Mach 
disc  to  the  contact). 

Total  number  of  grid  points  (including  points  on  both 
sides  of  both  shocks  and  the  contact  except  for  the  free 
stream  side  of  the  blast  wave) 

Axial  distance  (in  feet). 

Natural  logarithm  of  pressure  (in  psf). 

Gas  velocity  (In  fps) . 

Entropy  (In  ft^sec*-°R) 

Speed  of  sound  (In  fps). 

Gas  density  (in  slugs/ft^). 

Velocity  of  surface  of  discontinuity. 

Ratio  of  specific  heats  for  regions  1 and  2,  respectively. 
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FORTRAN  NAME 
CPI,  CP2 

RJET 

XHE 

PE 

AE 

PINF 

UINF 

AINF 

RHINF 

TIME 

TIMEF 

KK 

KO 

LL 

I PUNCH 


DESCRIPTION 

Specific  heat  at  constant  pressure  for  regions  1 and  2, 
respectively. 

Radius  of  gun  bore. 

Muzzle  exit  Mach  number. 

Muzzle  exit  pressure  (in  psf). 

Sound  speed  at  muzzle  exit  (in  fps) . 

Ambient  pressure  (in  psf). 

Ambient  gas  velocity  (in  fps)  (usually  0.0). 

Speed  of  sound  in  ambient  gas  (in  fps). 

Ambient  gas  density  (slugs/ft^). 

Current  time  (seconds). 

Final  time  (seconds). 

Maximum  number  of  time  steps. 

Current  number  of  completed  time  steps. 

Interval  in  the  number  of  time  steps  between  printing 
of  complete  flow  field  solution. 

Index  for  option  to  write  final  solution  on  TAPE12. 
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F0RTRAN  NAME 
OT 

OZ 

Def Ini tion  of  other 
of  their  usage. 


DESCRIPTION 


Time  increment. 


Increment  in  axial  distance. 


program  variables  should  be  self-evident  from  the  context 
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SECTION  III 
DESCRIPTION  OF  INPUT 


A.  PUNCHED  CARD  FORMAT 
Card 


Number 

Format 

Columns 

Description 

1 

15 

1-5 

starting  step  of  run 

15 

6-10 

final  step  of  run 

15 

11-15 

print  Interval 

15 

16-20 

output  file  creation  Index  (0  - no  output  file 
created;  1 - output  file  created  on  TAPE12) 

15 

21-25 

restart  indicator  (0  - initial  program  submittal; 
1 - program  reads  Initial  step  from  TAPE10) 

15 

26-30 

option  for  moving  origin  of  coordinate  system 
(0  - fixed  region;  1 - moving  origin) 

15 

31 -AO 

Input  unit  index  - specifies  the  units  of  the 
Input  (0  - metric;  1 - english) 

15 

A1-A5 

output  unit  index  - determines  the  units  of  the 
output  (0  - metric;  1 - english;  2 - non-dimensiona 

15 

A6-50 

acoustic  wave  index  (0  - normal  program  operation; 
1 - only  acoustic  wave  calculation  performed) 

2 

E10.0 

1-10 

starting  time  of  run  (seconds) 

E10.0 

11-20 

final  time  of  run  (seconds) 

E10.0 

21-30 

multiplicative  factor  on  maximum  time  step  cal- 
culated from  stability  theory  (recommended  value 
Is  1.0) 
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Card 


Number 

Format 

Col umns 

Description 

2 

E10.0 

31-40 

multiplicative  factor  on  minimum  grid  spacing 
as  determined  by  initial  number  of  grid  points 
(recommended  value  is  <_  1.0) 

E10.0 

41-50 

multiplicative  factor  on  minimum  grid  spacing 
to  determine  maximum  grid  spacing  (recommended 
value  Is  >.  2.0) 

3 

E10.0 

1-10 

radius  of  gun  bore  (nr  or  ft) 

E10.0 

11-20 

coordinate  system  Index  (0  - planar;  1 - cylindri- 
cal ; 2 - spherical) 

E10.0 

21-30 

ratio  of  distance  from  muzzle  to  origin  of  co- 
ordinate system  to  the  distance  from  muzzle  to 
the  Mach  disc 

If  the  restart 
not  required. 

i ndlcator 

is  equal  to  one  on  card  number  1,  card  number  4 is 

4 

15 

1-5 

last  data  point  in  region  1 (shown  as  1 R 1 in 
Figure  1) 

15 

6-10 

last  data  point  in  region  2 (shown  as  IR2  in 
Figure  1) 

E10.0 

11-20 

specific  heat  ratio  in  region  1 

E10.0 

21-30 

Specific  heat  ratio  in  region  2 

E10.0 

31-40 

specific  heat  at  constant  pressure  in  region 

1 (m*/sec*  °K  or  ft*/sec*  °R) 
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0—  Discontinuities 
Interior  Points 

Point  where  Streamline  Originates 


Hach 

Disc 


/ 


Contact 


Stream! i ne 


Region  1 


' Region  2 
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Card 

Number 

Format 

Columns 

Description 

4 

E10.0 

41-50 

specific  heat  at  constant  pressure  in  region 
2 (m*/sec*  °K  or  ft*/sec*  °R) 

5 

E10.0 

1-10 

2 2 
ambient  pressure  (N/m  or  Ib/ft  ) 

E10.0 

11-20 

ambient  gas  velocity  (m/sec  or  ft/sec) 
(typi cal ty  0.0) 

E10.0 

21-30 

ambient  speed  of  sound  (m/sec  or  ft/sec) 

E10.0 

31-40 

blast  wave  Mach  number  at  which  acoustic 
wave  analysis  commences 

6 

15 

1-5 

INUM-number  of  points  in  table  of  muzzle  exii 
properties  (maximum  of  25  points) 

7 

There  are 
properties 

"INUM"  of 
• 

the  following  cards  describing  the  muzzle  exit 

7a 

E10.0 

1-10 

time  (seconds) 

E10.0 

11-20 

pressure  (N/m*  or  Ib/ft*) 

E10.0 

21-30 

speed  of  sound  (m/sec  or  ft/sec) 

E10.0 

31-40 

Mach  number 

8 

15 

1-5 

INUMP-number  of  points  in  table  describing 
the  plume  centerline  Mach  number  distribu- 
tion 
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umber  Format  Columns  Description 

9 The  following  cards  describe  the  piume  centerline  Mach  number  dis- 

tribution. A maximum  of  35  points  are  allowed.  Input  the  location 
of  the  point  (ZZ)  followed  by  the  centerline  Mach  number  (ZMCL)  for 
that  point.  Input  4 pairs  to  a card  and  continue  on  the  next  card 
until  all  "INUMP"  points  have  been  input.  A sample  of  the  first 
ca;d  is  shown  below: 


9a 

E10.0 

1-10 

ZZ ( 1 ) - axial  distance  from  muzzle  in  bore 
radi  i 

E10.0 

11-20 

ZMCL(1)  - plume  centerline  Mach  number  at 
distance  ZZ ( 1 ) from  nozzle 

E10.0 

21-30 

ZZ  (2 ) 

E10.0 

31-40 

ZMCL(2) 

E10.0 

41-50 

ZZ  (3) 

E10.0 

51-60 

ZMCL (3) 

E10.0 

61-70 

ZZ(4) 

E10.0 

71-80 

ZMCl(4) 

If  the  acoustic  wave 
requi red. 

index  is 

equal  to  zero  on  card  1,  card  number  10  is  not 

10 

E10.0 

1-10 

the  location  of  the  contact  surface  (m  or  ft) 

E10.0 

11-20 

the  location  of  the  blast  wave  (m  or  ft) 

E10.0 

21-30 

the  initial  time  step  to  be  used  by  the  Runge 

Kutta  integration  routine  (usually  taken  to  be 
the  final  time  step  obtained  by  the  stability 
requirements  of  the  finite  difference  portion 
of  the  program)  (seconds) 
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Card 

Number  Format  Columns 


10  E10.0  31-40 


E10.0  41-50 


E10.0  51-60 


Description 

the  constant  in  the  Mach  disc  equation 
printed  out  at  the  last  finite  difference 
step 

the  constant  in  the  blast  wave  equation 
printed  out  at  the  last  finite  difference 
step 

the  constant  In  the  velocity  curve  fit 
printed  out  at  the  last  finite  difference 
step 


B.  Restart  Capability  - The  program  possesses  the  capability  of 
using  a solution  previously  stored  on  disk  or  tape  to  initialize  the  neces- 
sary variables  at  each  grid  point.  To  create  such  a file,  the  output  file 
creation  index  (card  1)  must  be  set  equal  to  one  (1).  When  the  final  time 
step  is  reached  the  necessary  variables  are  output  on  TAPE12  and  may  be 
saved  by  the  appropriate  use  of  control  cards. 

To  make  use  of  this  file,  the  restart  indicator  (card  1)  should  be  set  equal 
to  one  (1).  With  the  exception  of  card  2 which  must  be  changed  to  correspond  • 
to  the  time  step  saved  on  the  file  and  card  4 which  must  be  omitted,  the  pro- 
gram Input  deck  is  the  same.  The  variables  are  then  read  from  TAPE10  and  are 
used  for  the  Initiation  of  the  program.  At  the  completion  of  this  run,  the 
variables  can  again  be  output  on  TAPE12  if  desired.  Care  must  be  taken  to 
assure  that  the  files  and  tapes  correspond  to  the  desired  input  and  output 
data.  Typical  control  cards  that  illustrate  the  use  of  the  restart  capability 
are  shown  In  Figure  (2). 


CREATION  OF  OUTPUT  FILE 


JOB  CARD 
CHARGE  CARD 
REQUEST (TAPE  1 2 , *PF) 

ATTACH (LGO, DAWN, CY«1 , ID-ATLXXX) 

MAP (PART) 

LGO. 

REWIND (TAPE  12) 

CATALOG (TAPE 1 2, NAMESTEPXXX.CY-1 , ID-ATLXXX) 
EOR 


RESTART  USING  OUTPUT  FILE 


JOB  CARD 
CHARGE  CARD 
REQUEST (TAPE 1 2, *PF) 

ATTACH  ( LGO , DAWN , CY« 1 , 1 D-ATLXXX) 

MAP (PART) 

ATTACH  (TAPE  10,NAMESTEPXXX,CY=1 , ID-ATLXXX) 

REWIND (TAPE 10) 

LGO. 

REWIND (TAPE  12) 

CATALOG (TAPE  1 2 , NAHESTEPXXXX ,CY- 1 , 1 D-ATLXXX) 

EOR 

FIGURE  2.  CDC  SCOPE  2.1  OPERATING  SYSTEM  CONTROL  CARDS  ILLUSTRATING  THE  USE 
OF  THE  PROGRAM  RESTART  CAPABILITY 


TM  184 
SECTION  IV 

DESCRIPTION  OF  OUTPUT 


A.  Output  Format  - The  first  page  of  the  program  output  prints  a 
narrative  which  Informs  the  user  of  various  Input  parameters  such  as  the 
initial  and  final  time  steps,  print  interval,  free  stream  properties  etc. 

On  the  second  page  of  output,  the  input  values  of  the  muzzle  exit  conditions 
and  the  plume  centerline  Mach  number  distribution  are  printed  out  in  tabular 
form.  The  units  of  the  input  variables  are  specified  by  the  input  unit  in- 
dex as  described  previously  in  Section  III. 

On  the  following  pages  of  the  printout,  starting  with  the  initial  step,  the 
fluid  dynamic  properties  are  printed  at  each  grid  point.  As  a convenience 
to  the  user  the  locations  of  the  Mach  disc,  contact  surface  and  blast  wave 
are  specified  in  the  output  by  the  letters  MD,  C and  BW  respectively.  At 
intermediate  time  steps  where  a complete  output  is  not  desired  the  velocity 
and  position  of  the  discontinuities  are  printed.  The  program  has  the  capability 
of  printing  out  the  fluid  variables  in  either  metric  or  English  units  or  in  a 
non-dimensional  form.  The  output  unit  index  described  in  Section  III  controls 
which  unit  option  is  used.  A table  illustrating  the  unit  option  is  contained 
in  the  following  sub-section. 

The  form  of  the  output  for  the  acoustic  wave  computation  is  similar  to  that  of 
the  finite  difference  portion  of  the  program  except  that  only  non-dimensional 
output  is  used.  When  the  program  reaches  the  time  for  switching  to  the  acoustic 
analysis  (i.e.,  when  W(3)/a  < 1.10),  the  final  finite  difference  time  step  is 

output  along  with  the  necessary  parameters  for  resubmitting  an  acoustic  wave 
run. 
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B.  Identification  of  Variables 

STEP  - total  number  of  steps  taken 
TIME  - total  elapsed  time  (seconds) 


Output  Unit  index 

Description 

0 

1 

2 

z 

axial  distance 

meter 

feet 

Z/D 

p 

pressure 

N/m2 

lb/ft2 

P/P 

CO 

u 

veloci ty 

m/sec 

ft/sec 

U/a„ 

00 

s 

entropy 

2,  2 o„ 

m /sec  K 

ft2/sec2  °R 

(S-Sr«f)/Cv 

A 

speed  of  sound 

m/sec 

ft/sec 

00 

RH 

density 

kg/m^ 

slug/ft^ 

p/pm 

TEMP 

temperature 

°K 

°R 

T/T, 

MACH 

Mach  number 

- 

- 

- 

W(1) 

velocity  of  Mach  disc 

m/sec 

ft/sec 

WD/a^ 

W(2) 

velocity  of  contact 

m/sec 

ft/sec 

W(2 )/a 

oo 

V(3) 

velocity  of  blast  wave 

m/sec 

ft/sec 

w(3)/a 

OO 

2(1) 

location  of  Mach  disc 

meter 

feet 

Z(l)/D 

Z(2) 

location  of  contact 

meter 

feet 

Z(2)/D 

Z(3) 

location  of  blast  wave 

meter 

feet 

Z(3)/D 

C.  Sample  Input  and  Output 
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TYPICAL  INPUT  IN  CARD  IMAGE  FORMAT 


TM  184 


MUZZLE  BLAST  ANALYSIS 

STARTING  STEP  * 0 - - 

— -FINAL  STEP  s JO  

PRINT  INTERVAL  =10  - 

STARTING  TIME  s 0.  SEC,  --  - 

FINAL  TIME  = .30000E-02  SEC,  - 

LAST  PT.  REGION  NO,  1 = 7 - --  - 

LAST  PT.  REGION  NO.  2 = 13 

MOVING  COORDINATE  SYSTEM  OPTION  = 0-  

RADIUS  OF  THE  JET  = .91300E-02  FT,  

FLO*  INOEX  * 2 

- — ORIGIN  OF  SPHERICAL  SYSTEM  s .400E*00  TIMES  RJET  

FREl  STREAM  PRESSURE  * ,2120QEfQfl  LBS/FT/FT  - — 

FREE  STREAM  VELOCITY  = 0,  FT/SEC  - 

.FREE  STREAM  SPEED  OF  SOUND  « .UZOOEtOa  FT/SEC 

SPECIFIC  HEAT  RATIO  (REGION  1)  =1,25  

SPECIFIC  HEAT  RATIO  (REGION  2)  = 1,40  - - - 

SPECIFIC  HEAT  AT  CONSTANT  PRESSURb(REGION  1)  * 9507.1  FT**2/SEC* *2-0EG.R 

SPECIFIC  HEAT  AT  CONSTANT  PRESSURE  ;‘>FGION  2)  c 6006.0  FT **2/SEC**2-DEG.R 

BARREL  EXIT  PRESSURE  = .l29SlEt07  LBS/FT/FT  

BARREL  EXIT  MACH  NO.  = ,10000E*01  ...  — 

BARREL  EXIT  SPEED  OF  SOUND  « ,42«93E*0«  FT/SEC 


FIGURE  A.  FIRST  PAGE  OF  PROGRAM  OUTPUT  - 


i 


-18- 


TH  18* 


i 


i 

I 


ac 


x 

u 


Ml 

JC 


' tr 

l Ml 
I ►— 


u 


I 

I 


tr 

•v 

X 


ooooooooooooooooooooooooo 
♦ ♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦ 
||J  m Ml  Ui  Ml  Mi  MJ  Ml  Ul  Ml  j|  ||J  dj  Ml  Mi  Ml  Mi  Ua  Ml  Ul  Ml  Ui  IaJ 
O N 90NM\|NlMM13SeiONOO  OOOOOOO 

o^jim®  a «n  t'j  O'  >o»  ,'i-«NrO'<wi^rv(\iMOBOo 


OOOOOOOOOOOO 
♦ ♦♦♦♦♦♦♦♦♦♦♦ 
Mi  mJ  Ml  Ml  mI  Mi  Ui  mI  Mi  Ml  ^ 
SrM99909-0«<fl«e 

O OxO  3 S 


II  : i 

M«4>«rMMr\iiM'MrvjrgiMru 

OOOOOOOOOOOO 

nj  ^ m ^ 

O <0  9 000000000 

J"»KlinOOOOOOOOO 

3NB  O - (\l  3 DCOOB 

BNC^x-x--ixMOB 


f 


Ml 

X 

3 

-I 

a. 


t 

1 


I 

i 

» ! 

i 

1 

1 


1 


000000000000000000000 

J ^ ^ yj  ^ jj  jJ  ^ ^ ^ yj  |^j 

000  o 0000000000000c  0000 

9000000000000000000000 
X OOOOOOOOOOOOOOOOOOOOO 


9 9 9 9 9 9 9 
OOOOOOO 
♦ ♦♦♦♦♦•♦ 
(4)  ^ ^ 14J 

3 lO  x l\i  111  O iT 
ry  » N X (M  o 
9 K*  t*\  r*>  PO  9 » ?n 


9 9 9 9 9 9 9 
OOOOOOO 
♦ ♦♦♦♦♦♦ 
yj  I4]  (4J  yj 

N033ffCd 

O ® D r*-  r>»  r>*  f>. 

!\i  m f\i  r\i  r\i  oj  r\j 


9999999 

OOOOOOO 
♦ ♦♦♦♦♦♦ 
^ ^ ^ ^ llj 

r»  iriftrtio 

oss 

MI4««^0 

(M  ryj  rvi  fvj  ru  <M  r\j 


*■»  ^ f«-  f«»  <0  -C 
000000 
♦ ♦♦♦♦■♦ 
Ml  Mi  .Mi  Ml  Ml  Mi 
00  B 'V  N N X 
O N.  9-  «*  •*  9 

im  — o © in  o 

x*  x x ^ B 


OOOOOOOO 
♦ ♦♦♦♦■♦♦♦ 
Ml  Ml  Ul  Ml  Ml  Ml  Mi  U 

r*.  -n  *n  9 -•  y s* 

MD039lHax 

9»9^fy^»n^ 

B ^ N < D IfllT  9 


B « B B Bl/un 
OOOOOOO 
♦ ♦♦♦♦♦  + 
yj  14J  l4l  ,jj  yj 

^ ^ X x «x  N 

in  *n  in  » in  c m 
BIU»B9(MM 
3 3 Ifl  i'l'O  B »0 


9 

OOOOOOOOOOOO  90900000 

1111111111111(111111 

yj  |^j  mJ  4J  MJ  Ml  ij  mJ  Ml  !ii  Ml  Ml  ^|  mJ  i|  Mi  Ml  Mi  Ml  lii 

X OOOOOOOOOOOOOOOOOOOO 
•4  009f\|x^NB^x»KB7MOOBOO 

^ ®B»999®,n®>n»ru^m»*.fM^.-4M)oo 

9>xxM!\|414|  3 J^in  a BNNBB'M'O 


1 


1 


18 


FIGURE  5.  SECOND  PAGE  OF  PROGRAM 
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FIGURE  6.  TYPICAL  OUTPUT  OF  FINITE  DIFFERENCE  CALCULATION 
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FIGURE  7.  TYPICAL  OUTPUT  OF  ACOUSTIC  WAVE  ANALYSIS 


TM  1M 
SECTION  V 
PROGRAM  LISTING 


PROGRAM  DA*NA( INPUT, OUTPUT, PUNCH, TAPfc5«lNPUT,TAPE6«0UTPUT, 

JJTM5E10j_TAPE12)  

COMMON/ A/2 (200), P(200)#U(200),RH(200),S(200)# A (200) »*(3) 
COMMQN/B/Z N ( 2 0 OJ ,JPN_(2 Q0> i UN (200)  §RHN(200) »SN(200)«AN(200)»*N(3) 
COMMON/ 0/ 1R(2),GAM,xJ,CP,0Z,0T 

CQMMQN/F/RJ&T, XME,PE, AE,PINF, AINf  .UINF.RHINF 

COMMON/G/KK,LL,KO#TIM£,TIMEF 

CQMMON/IPu/  1PUNCH 

COMMQN/ZZERO/ZZtRO, !MOvt,OZZERO#OTZZER,PACT 

COMMON/G_A/GAMi,GAM2iCPl#CP2  

COMMON/BLAST/ IB*, b*CON, SPEED, RJET2 

m^PN/.DSHUU/lOLBuG 

COMMON/UNl TS/IuNI 7 , 1uDUM,1uNQUT»IUOUT»FTME,PuNIT,DEUH*,PhUNIT 

Common/ AMI /3. #maC  * 

OATA  H1/2M  /#h2/2hMD/,HJ/2H  C/,HA/2hB*/ 

PI*CL*fi 

CALL  INDATA 

FAC1NF«CP2«{GAM2M.0)/CAM2 

FACT2«FACINF 

FAC_LUCPi*J  GAMIM.OI/GAMI 

3 CONTINUE 

IHlSlRtll 

IR|P*IR1 ♦ 1 

\*i*l*U) 

IP2P«IR2tl 

IF(KO.GE.KK)  go  TO  a 

IF((KO/LL)*LL.NE.KO)  go  TO  5 

A CONTINUE 

NRITE(6,7)  KO#  TIME 

T FQRMATtlHl, 10X,*STEP  »«, IA, 10X,«TIME  »»,E 1 3.5, 1 X , AHSEC. . // ) 

PR1TC (6, 8 ) 

8 FORMA T ( A * * I * 1 1_X  ,_*Z *,1  3 X , *P » , 1 iX , , 1 3X , »S* , 1 3X  , « A « , 1 ?X , «RH« , 

U IX,* TEMP*,  iOX,*MACH*) 

IULIP(2)*A 

FOUM«AINF*TIME 

DO  9 I« 1 , II 

FACRtGsFACTl 

I Fd.GT.IR  1 J FAC  RfcGsFAC  T2_ 

H5*H1 

IFCI.EQ.l .OR.I.EQ.2)  H5*H2 

I F ( I , EU . I R 1. OR . I. EU . I R 1 P ) H5*M3 

IFC1.GE.IR2)  H5«HA 

GO  TO(22,2A,26), 1UOUT 

22_C0NTJNUE  

25*2(1 )*FTME 

P5*EXP(P(I))*PUNIT 

U5*U( I ) *F TME 

A5«A ( I ) *F TME 

R5*RH( I ) aRHUNIT 

85*8_(I)*FJME**2tqEGRI( 

TEMP«(P5/(R5*FACMEG))/FTME**2/DEGRK 
GO  TO  50 
2A  CONTINUE 
25*2(1) 

P5*EXP(P ( I ) ) 

U5*U( I ) 

A5«A(i) 


23 


TH  184 


RSsRH(I) 

asi&rn 

TEWR*(P5/(R5*F  ACREG) ) 

CO  TO  SO 

26  CONTINUE 

2S»2<I1/RJET2 

P5*EXP(P( I ) J /PINE 

U5«Utn/0INF 

A5«A(I)/AINF 

• RS«RHf IWRHINF  

TEMP«(PS/Rb) *F AC  INF /F ACREG 

l£U.6I.lRn_6U-ICLifl 

S5«GAM1*(S(I)-SU))/CP1 

60-IP  ..40 

30  CONTINUE 

S5«GAM2*(S(I)- SU«  2PD/CP2 

40  CONTINUE 

50  CONTINUE 

XM5sU(I)/A(l) 

ISP-1IE(6,10)  liH.5#25j.P5LtUS<A&i»5jR5iI£MPjXMS 

2000  FORMAT (8X#5E14,5) 

j-COMlPUE 

10  FORMAT (I5«1XA2#BE14(S) 

mR1TE(6» 12) 

12  FORMAT!/) 

lElKfl«6E«gM-6fl-TQ...fc 

R«ITE(6,21) 

2 1 F QpMAJliiJL*JlUJ.« IJQX t*12 J * ALOX* Z111*1PA*112 ) *i0x*£L3)*  - 

llOX*TIME*,llX,*DTt) 

5 CONTINUE 

21*2( 1 J/RJET2 

22«2(IR1)/RJ£T2 

23*2(IR2)/RJtT2 

IF(Z3.GT.l .0) 

1BNC0N«23  *(1.0*(AINF/m(3))*«2)*SQRT(ALOG(23  )) 

IF(23.GT.  1 .0)  BwTCON«TIMEtl.0E»Q6»CRJET2/AlNF).*tBwCJN» 

15ORT(AlOG(2i))-23)*l.0EtO6 

2MP«Q.69«iXME  * S U H L(  G A Ml*P£/£  XJM  P LtU  1 

GO  TO(32,3«.3b), IUOUT 

32  CONTINUE 

WlOUMSN(l)*FTM£ 

N2PUM«»(gj»pTME 

W30UM*b(3)*FT-Mt 

210UM«2(1)»FTME 

220UM«2(IR1)*FTME 

230UM«2(IR2)*FTME 

GO  TO  39 

34  CONTINUE 

W1DUM>*(I) 

W20UM«w(2) 

W30UM«N(3) 

21PUM«2(1) 

22DUM«2(1H1) 

Ilmiilliigi 

GO  TU  38 

IVCONUNyg 

610UM«N(1)/AINF 


ft2DUMs*C2)/AINF 

w3DUMkh( 3 ) / A 1 N F 

Z1DUM*ZC1)/RJET2 

Z20UN*ZCIR1)/RJ£ T 2 

23DUM«Z(IR2)/RJET2 

36  CONTINUE 

UNITE (6, 1 1)  *tOUM#w2DU*#*3DUM,ZlDU*rZ2DU*'Z3DUM,TIME'DT 

11  FORMAT ( S X »J)E 1 4 . 5 , 2F8.2) 

IFUBH.EQ.O)  60  TO  la 

_ LL5LL0UM  

ZCONcZ(INl) 

ZBwsZ  ( 1W2 1 , 

RCT*Z ( INI P)»FOUM 
SGNsSIGNC 1 1 0»  RCT ) 

00  60  1*1  HIP# IR2 

__RCTsZ(l)-FOyM 

SGN1*3IGNC1.0,RCT) 

1FCSGN.EQ.SGM)  GO  TO  SR 

IRCT«1 

IRCTHkI-1 

RCTM«Z(IRCTMJ-FOUM 

TERM* (Z( I R C T ) / N JE  T 2 ) * *J  * .CRH (IRC  T)/RHINF»1 ,Q«UC 1RCT ) /SP EEO) 

TERMM«(Z(IRCTM)/RJfcT2)**2*(RM(iRCTM)/RHlNF-l,Q*U(INCTM)/SPEE0) 

RATRCT«RCT/CHCT-RCTH) 

RMORCT*TEHMfRATRCT*( TERMM-TERM) 

GO  TO  70 

59  CONTINUE 

SCN*36N1 

60  CONTINUE 

RHQRCT«1.0 

70  CONTINUE 

WRITE (6,1 000)  CQNMQ, BwCON, RHORC T 

1000  FORMAT!//,  IX, *TH~E  CONSTANT  IN  THE  MACH  DISC  EQUATION  IS*,E13.5, 

1//,1X,«T HE  C ON S T A N T IN  THE _BL AST  NAVE  E QUATION  IS*, El 3 . 3 L 

1//, 1X,*THE  CONSTANT  IN  THE  VELOCITY  CUHVE  FIT  I3*,E13,5) 

CALL  ACOUS(ZCON,ZBw,DT,CONMQ,RHORCT) , 

14  CONTINUE 
CALL  SETN 
CALL  STEP 

_ ___CALL__5H0CK  ( 1 ) 

CALL  SHOCK (3) 

CALL  CONTCT(2) 

CALL  INT  PT 

IF(IMOVE.EQ.O)  GO  TO  20 

ZZER0*ZN(1)*FACT 

DZZERUswN(l)*FACT_ 

DTZZER«(ir,NU)-w(l))/Of 


20  CONTINUE 


Il*IN(2)fl 
00  1 1*1,11 

Z (I)*Z  NCI) 

P <l)*P  NCI) 

U CI)*U  N(I) 

3 C!)*S  NCI) 

A C I ) *A  N ( I ) 


1_RH(  I ) *RHN  ( I ) 
00  2 1*1,3 


TlM£*TIMfc4DT 

If  ( ( * ( 3 ) / * 1 NF ) . Ulu fljfMEM J-.1ML1 

IF(IBW.EQ.O)  CO  TU  5 

CALL  MUZZLE (T  IME*PEXIT*AEXIT«XMEXIT) 

CONPDsZ(l)/(XMExIT*SQRT(CAMl*PEXlT/PINK))/KJtT2 


LL«i 

G0_LQ._I  _ 

6 CONTINUE 


(a*  1 JUiMI 


TH  184 


subrout i n t cluti.bt) 

COMMON/ A/Z ( 200 ),P( 200) .UC200) .RH(200).S(200). A(200).w(3) 

COMMUN/B/ZN(200) ,PN(2O0) ,UN(200) ,RHN(200),SN(2QO) , AN(200), **(3) 

C0MP0N/D/IR(2)jUAM,xJ,CP,0ZjDT 

COMMON/F/RJE1 ,XNE,PE# AE.PINF, AlNF,UlNF,KHlNF 

COMMQN/GA/  GAM  1 , b,» M2 , C P 1 . CP2 

COMMON/G/KK,LL#KU, TIMfc, TlMtF 

C0H10N/wCL/ZZ(3S),Z'1CL(35)»EMCL(3SJ,INUMP 

COPm6n/IN1T/PM(25),am(2S)»XMM(25),TM(25) 

CJjMMON/DEBuG/IDEbUG  

DIMENSION  Z«U)'XM(1),DXM(1),02X«U) 

OATA  MAX. ERRVAL/1 . 0.001/ 

ZEsO.O 

I T E R * 0 

RAT  10*0 , 5 

litfilj  10.0 .0 ) _ UN  t U*U.UJ 

T INI T*T IME+DT 

HUJlLl  3 

IF(IDEBUG.EU.b)  *RITE(6.1)  I T 1 » ZE , RAT  10. T IME , DT , T INI T * T1 , BT 

<LQl  COM-lNUt  

TDUMSsRAT IO*UNlT+(  1.0-HAT  1 0)*T1 

C ALt_MUZZLE  ( TUU.*iS  , PE , AEj  XME) 

UE*0.5*(XME*AEtUN(l)) 

TDUMP8T1NIT-(ZN(1)-ZE)/UE 

ERR=(TDUMS-TOUMP)/TDUHS 

JMlOEBUG.EQ.b)  Nfi  I TE  (6r  1 ) I TER , TDUMS,  TDUMP, PE.,  At , XME , UE  f ZN ( 1), ERR 

1F(TDUMP.LT,U.0H,  IDUMP.GT.T1NIT)  GO  TO  100 

JF(ABS<ERR)tLTtEHWVAL)_GO_TO  300  

ITEH*ITER+1 

1F(ITER.GT.20)  GO  TO  200 

RATIO*  (1 1-TDU«P)/(T1-TIMT) 

GO  TO  20 

100  CONTINUE 

*RITE<6,  1000)  

1000  FORMAT(lXf*STREAMLlNE  IN  SUBROUTINE  CL  IS  OUT  OF  BOUNDS*) 

CALL  EXIT 

200  CONTINUE 

NRITt(6,2000) 

2000  FORMAT ( 1 X # *TUO  MANY  ITERATIONS  IN  SUBROUTINE  CL*) 

CALL  EXIT  

300  CONTINUE 

GAMgGAMj  

1*1 

ZZl*ZMl)/RJETMt  

GlsGAP-1, 

XME2sXMfe.xME 
Fl*l . 4Gl*XME2/2, 

GisGAM/Gl 

G«*l ,/Gl 

ZM(1  )*ZZ1 

CALL  SPLINT(ZZ»ZMCL, INUMP»ZM»MAX,xm.DXM.D2XM,£mCL) 

XMCLsXM(l) 

XMCL2*XMCL*xmCL 
F2si .♦Gl*XMCL2/2. 

PPE«CFl/F2)**Gi 
RHRE*CF1/F2)**G«_ 

PN(I)sPPE*Pt 


27- 


ItULSi 


PHf=GAM*PE/At/AE 

RHN(I)sfiHRE*RHE 

AN(I)sSURT(GAM*PN(l)/RHN(D) 

UN(1)SXMCL*AMII)  

SN(I)=S(l)fCPl*((ALOS(PN(in-P('n)/GAMi-A'LUG(RHN(n/RH(l))) 

PN  ( I ) = ALUG(PN (1)) 

IF  CIDfcBUG.Ed.~t>)  wRITE(b#l)  I # PN  ( I ) # RHM  ( I ) , AN  ( I) , UN(  I ) , SN  ( 1 ) 

1S(I)#P(1)#RH(1)  

RETURN 

1 FCRK  AT ( 1 I b * Rfc 1 i.b) 


TM  184 


SUBROUTINE  CONTCT(A) 

__CP«MON/A/Z(?00),R(200),U1200),RH(200),S(200), A(200),*{3)  _ 

CO~MMON/B>2N(200),PN(200),UN(206')'RHN(200),SN(200)'AN(200)i*N(3) 

C OMMON / C ✓ Z I < 2 0 0 ) , P I ( 2 0 0 ) * U I ( 2 0 0 J »RHI(200),SI( 2Q0),A1(2 0 0 j 

COMMON/D/IR(2) ,GAM,xJ,CP,PZ,PT 

cqmmqn/e/bt,al 

COMMON/GA/  GAM1,GA^2,CP1,CP2 

COPMON/ZZERO/ZZERU, IMQvE»OZZERO»OTZZER»FACT 

1 • 1 R Cl ) 

I T | = I 

AL*1. 

BT*0. 

MN(K)*W(K) 

_ _UN(I)*»VN(K) 

uN(itn*UN(n* 

2 ZN  ( I )*  Z(I  )t  ( AL*  *{KUBT«wN(K))*PT 

i ZNiiti)*ZN(ij 

UlCNN(K) 

U2*U1 

GAPcGAMl 

CP»CP1 

CALL  LPUI_NT(I,  

GAMCGAM2 

CP*CP2 

CALL  LPUINT(H1>*1«#*1#2) 

_ HN*ZN(1)-ZZERU _ __ 

Rl*ZI(l)-ZZtRO 

R2*Z1(2J-ZZEH0  

A1*GAMI*(AL/A1C1)%BT/AN(D) 

B1«GAM2«CAL/AI(2)»B1/AN(I»1  n 

IFdPOVE.LU.O)  GO  TO  200 

_ A<I«GAM1QJ»(ALo(UIU)»DZZEHO  )/*t»BT*IUN(l)»PZZKKO  )/«N)_ 
A4«A4*PT 

BA«GAM2»X J*  C AL» (UI (2)»PZZERQ  >/R2»BT» (UN( I»j )»PZZERQ )/RN) 

B«cUa*DT 

GO  TO  210 

200  A4*(AL*Ul (1)/H1*BT«UN(I)/KN)*0T*XJ*GAM1 
B4s(AL*Ul (2)/H2  + BT*UN(IM  >/RN) *X J*GAM2*PT 
1T0  PCl*PI(l)-lUi-Ul(l))*A|-A« 

PC2«Pl(2)*(u2*Ul(2))*6t-B4 

Apcmxptpcn 

APC2sfcXP(PC2) 

ER*(APCl»APC2J/(APCltAPC2)*2.0 
IF(ABS(EH).LT,l.t-0«)  GO  TO  7 

in*  i mi 

IFdTl.Lt.  IS)  GO  TO  777 
WRITE (6# 111) 


111  FORMAT  (5X » »TUU  MANY  ITERATIONS  I'M  SUBRQUT  iNt  CUNTCTo) 


» 

STOP 

777  IFdTl.GT.2)  GU  TO  ltt 

BET«WN(2) 

ERlsER 

BETlsBET 

BETsl  .01 *BET 

GU  TU  15 


14  PURisBET 1»ER1*(BET»BET1)/_(ER»ER1) 
ERlsER 

- . -29- 


L 


BETlsBEI 

BETsDUM 

15  *N(K)sBET 
ZM 11 5111 

)»aN(K))*PT 

GO  TO  1 
N 


PN ( I ♦ 1 > *PN( 1 ) 

-UMUfUi 

UN(Ifl)*U2 

J&Ll  )•*<!)  

sn(i  + i )=s(  i+n 

kh(t  n.p(i  ) /gan 


AN(l*l)*SQKr(GAM2*tXP(PN(lMn/RH*(I*l  ) ) 

If(BT.GT.O.)  GO  Jfi  9 

BT*,5 


RHN(I  )=PN(I  )/GANl-(SN(I  J-S(I  ))/CPl*CK 

_RHNli)=£XP(RHN(I)) 

CKsALOG(RH(I  + l ))«P(U1)/GAM2 

RHNUM)=PNUM)/GAM2-($NU4n*SU.fU)/CP2*CK 
RMN(1*1 )=EXP(RHN(1+1 ) ) 

AMUsSURTIGAM 


9 RETURN 
END 


TM  184 


SUBROUTINE  FS 

CQNMON/B  /2  N ( 2 0 OJLi  PN (200 ) « U NJ  2 00),RHN(2  Q 0 ) , SM  20  0 ) , A N ( 20  (Uj  * N(3) 
COMMON/O/ 1« ( 2 ) # GAM, X J , CP , DZ  # 0 T 

COMMGN/F/RJET jXME,PEi AE.PIN F.  AlNF  . UINF . RHINF 

COMMON/GA/  GAM1,GAM2,CP1#CP2 

UlRilJ  t.l 

PN(I)rPINF 

UN(I)sUl N F 

RHN(I)sRHINF 

_*n_iu*o, 

AN(I)*SG«T(GAM2*PN(IJ/RMN(I)) 

pn(I)«alugcpn(1)) 

RETURN 

END 


- 4 


nnnl 


TMI8J* 


SUBROUTINE  INDATA 

COMMON/ A/ Z (200),P(200),U(200),RH(200)»S(200),A(200)#rt(i) 
COMMON/B/ZN(200) ,PN(200),UN(200),RhN(20Q) ,SN(20Q),An(2u0),  *N(3) 

C0MM0N/D/IR(2),GAM,XJ,CP,DZ,DT  _ _ 

COMMON/F/RJET,XME,PE,AE,PINF, AINF,UINF,RHINF 

CQMMQN/G/KK , Li. » * 0, TI HE , T IMEF 

common/dz/ozmin,ozmax 

COMMON/ DT  STAt*/0TSTAS  _ _ 

COMMON/ IHU/  IPUNCH 

COMMON/GA/  GAV1 , GAM2,CP1 ,CP2  

COMMON/ F 1 1 /EMP (25 ) # EMA (2S) » tMM( 25) # INUM 

COMMON/ 1 N1 T/PV ( 2b ) , AM ( 2b) , XMM( 25 ) . TM(25)  

COMMON/MCl/ZZ(ib) ,ZMCLl35)#tMCL(3S) , I NUMP 

COMMON/ Z ZERO/ Z ZERO, I MOVE, DZZERO# DT ZZER,  F AC T 

COMMON/RRAT/PRAT 

COMMON/OEBUG/1UEBUG 

COMMON/BLAST / 1 Ba , Bi/vCON,  SPEED,  R JET2 

COMMON/BwANAL/BaMACH 

COMMON/UNITS/IUNIT,  I JOUM,  IUN'JUT,  IUUUT,FTMfc,PUNlT,i)EGR^,RHU  *11 

DIMENSION  T 0(1), PO(l), DPO ( l ) , 02PUC  1 ) 

DIMENSION  M 3 ( 2 ) # m u ( 2 ) , M5 ( 2 ) # H6 ( 2 ) , H 7 ( 2 ) 

OATA  H1/«HRJET/,m2/«HZ(1)/ 

DATA  I DEBUG/ 0/ 

DATA  FACP,FACA/1. 0,1.0/ 

DATA  FTMt, PUNl T,DE&RK,RMUNIT/0.30a8,« 7. BB02b8, 1.0,515,379/ 

DAT  AH8/4HSEC , / 

DAT  A H3(  1 ) ,•  HU(  1 ) , H5(  1 ) ,H6(  1)  »H7  ( 1 j / JHM'  ,9hN/M**2  ,6HM/SEC  V 
U0HM**2/SEC»*,8H2-DEG.K  / 

DATA  H3(2),Ha(2),Hb(2),M6(2),H7(2)/3HFT , , 9HIBS/FT/F T , b*F T /SEC 
1, 10HFT»*2/SEC»,HH»2-DEG.R/ 


READ  INPUTS 


_ READ (5# 1 ) KO,KR,Lt/ IPUNCH, IRSTRTriMOvE, lUNIT^IUNOUT, IB* 
lUDUMclUM  T ♦ ! ‘ 

lUOUTslUNOUT+1 

HDUMcHl 

IF ( IMOVE • EU, 1 ) HDUMch2 

READ (b, 2)  TIME, TlMEF,OTSTAd,DZMlN,DZMAX 
REA0(5,2)  RJET,XJ,FACT 
IFUUNIT.EQ.O)  RJETrRJfcT/FTME 

I F ( I RSTRT  « GT  « 0 ) GO  _T_0_8 

RE AD (5, 3)  IH(1),IR(2),GAM1,GAM2,CP1,CP2 

IF ( IUNI  T ,EU, 1 ) GO  TO  8 

CPi*CPl/F  TME**2/DEGRX 
CP2bCP2/FTME**2/DEGRK 
6 CONTINUE 

Rt  AD  ( b»  2 ) PINF,U1NF,AINF,BwMACH 

IMlUMl.EO.i)  GU  TO  11 

PINFsPlNF/PUNIT 

UINF*U1NF/FTME 

AINF*AINF/FTME  *_ 

' FACP*PUNIT 

FACAbFTME 

II  CONTINUE 

READ(b,  1 ) INUM 

DO  10  1*1,1 NUM 


HUB* 

READ (5, 2)  TM(I)#RMCI),AM(1),xmm<I) 

IFCI.EU.l)  TM u si m lij 

PM(I)sPM(I)/PUNIT 

APm«AMtn/FTME 

10  CONTINUE 

£A kL_5 PLlN EXT PM#  j \uMx£.«P  1 

CALL  SPLINE(TM,AM,INU*»EMA) 

CALLJSPLINE  ( 1 v’#  I 

CALL  MU22LE(TlMt,Pt,AE,XMfc) 

READ (5« 1 ) 1NUMP 

READ (5, 2)  (ZZ(I)»2MCL(I)#I*lf INUMP) 

JFtlBN.EQ.O)  GO  TO  1 

READ (5, 2)  2C0N,Zttw,DT,CONMD,BWCON,RMORCT 

1F.IXUN II » t Q,  1 ) GO  JU0_2 

ZBW«Zbw/F TMt 

2C0N8ZCUN/FTME 

7 CONTINUE 

IEJL1R  5I«I.  GI  

9 CONTINUE 

RHINF8GA^g»PlNF/AINF/AINF 

SPEEDsAINF 

RJE  T2=R jLT*2j 

IF(IBN.EQ.l)  call  ACUUS(ZC0N,ZBw,DT,C0nMD,RH0RCT) 

LFJJUNIL.EIM)  Gl‘  T^>  320 

CP1DUM*CP1»FTME**2*DEGRK 

iP2.DUMsCP2*FTMt**2*DEGRi< 

PEDUM=PE*PUNIT 

AEDUM»AE*FTM£ 

RJDUP*RJET*FTML 

PlDUPgPlNF  • p u n I t 

UIDUMsUINF*FT^t 

AlCLUMs  A J NF  * F TM  E 

GO  TO  325 

tea.  continue 

CPIDUMsCPI 

CP2DUN*CP2 

PEDUP*PE 

AEOUMsAE 

RJDuMspJET 

PIDUPsPIMf  

UIDUPSUINF 

AJOUMsAINF  _ 

125  CONTINUE 

C wNlTt  INPUTS 

WRITE (6» 1652) 

WRIT_E(6,200)  WO«KKfLL»TIMEiH8»TIMEF»Hd»IR»  I MOVE 

J*’XJ4.5 

__  WRITF.(6,20l)  RJOUM,  H3(1UDUM)»  J»FACT,HDUM 

WRITE  (6, 202)  P IOU*#  Hit  { IU0UM)  , UIOUM,  H5(  IU0UM)  , A IDU",  H5  ( IUD  J* ) 

WRITE  (6. 20b)  (,AMtfGAR2,CPlDUM«Hb(lUDUM),H7(IUOUM)  > CP23UM, Hb ( 1 UDUM) 

1,H7(IUDUM) 

_ WRITE  (bf  207)  PEDUNjJHA^ IUDUM)  t XmE # AEDUM*  H5( JUDUNJ 

WRI TE (b#  205) 


33 


o n 


TH  lflfr 


*R1TE(6,2Q4) 

INUMQsINUM  

IMINUMH,LT.  1NUMQ)  nUMUsINUMP 

D0  326  1*1#  I Nu**Q 

POUHMysPMt I)*FACP 

ADUMMY=AM ( I) «F  AC A 

WRITE (6# 500)  T*l(l),POUM*Y,AQUMMY,XMM{I),ZZ(l)#ZMCLCI) 
326CONTINUE 

INUMO*lNUWQtl 

IFUMJV-1MJMP)  307,308,309 

307  wHITt(6,501)  (ZZ(1)#ZMCL(1)#I=1^UMU#INJMP) 

GO  TO  30b 

309  CONTINUE 

DO  327  IsINUMU, INUM_  

PDUFfcY*PM(l)*FACP 

ADUMMYsAM(l)*F ACA  

WRITE(6#502)  T C I ) # PDUMMY#  AOUMMY# x~MM(  I ) ' 

327  CONTINUE 

308  CONTINUE 

DO  12  1 = 1, INUMP  

zzTi>*zz(mi.o 

12  CONTINUE  _ _ 


INITIALIZE  FLO*  PROPERTIES 


CALL  8HL1NE(ZZ#ZMCL,INUMP,EMCL)„  

IF(IBSTRT.GT.O)  GO  TO  3«2 

CALL  INI T 

1»1*IP(U 

1W1 1 sIRl ♦! 

IR2=1R<2) 

_ IH21*IR(2)+1  

IR«=IB1-1 

I R5* I R 1 ♦?  

IR6«IR2-1 

P2*EXP(P(2) ) 

PIR1=EXP(P(IR1)) 

CK=ALOG(RH(l))-P(l)/GAMi 
DO  330  I=3#IR« 

RAT  =FLOAT(I-2)/FLOAT(lRl-2) 
Z(I)sZ(2)t(ZClRn*Z(2))*HAT 

U(I)*U( 2)  ♦ ( UJJRJJ  - UJ2 ) ) *RAT 

S(I)sRAT*AL0G10(S( IR1 )/S(2  )) 

xn*io. _ _ 

S( 1 )=S(2  )*S(I) 

P(1)=P2+(PIR1-P2)*WAT 

RH(l)  = ALOGlPCI))/GAMl-(S(I)-S(l))/CPl*C*v 

Hh( I )=E  xP(Rh( I ) ) 

A(I)sS0RT(GAMl*Pll)/«H(I)) 

P(I)*ALOG(H( I ) ) 

330  CONTINUE 

P2*EXP(P( IRll ) J 
PIR2sEXP(P(lR2)) 

CK«AL0G(RH1NF)»ALUG(PINF)/GAMZ 

DU  340  I = 1H5  # IR6 

RAT  *FL0AT(1-IR11) /FLOAT (IR2-IR11) 

/tiy*Z( 1R1 l )♦ C Zl 1R2)»Z ( I W1 1 ) )*RAT 

.34. 


TM  18  k 


U(I)sU(IRU)*(U(IR2)-U(lRU))*RAT 

8U)sRA I_a  A L QG  1 0(S(IR2)/S(1R11 ) ) 

SCI)*10a**S( 1 ) 

111M1  H.L11  * 5 1 1 1 

P(I)sP2*(PIR2-P2)arat 

RHCl)sALQG(P(l))/bAMg-(S(l).S(IR2in/CP2»CK 

RH( I )s£xP(RH ( 1 ) ) 

ALLUSOM  TXG  A M2*£UJ /RH(l)j 

P ( I )*ALOG(P( 1 ) ) 

3RQ  CONlIf-Uk 

i«2  CONTINUE 

RETURN 

300  CONTINUE 

C ******************  ******** *************************** 
C READ  RESTART  VARIABLES 

READ(IO)  IK l 1 )»IR(2)f SAMI, GAM2, CPI  rCP2 

READ(IQ)  (*(N1,N«1,3) 

READ(IO)  ZZERO,OZZERO,DTZZER 

IiSULl21±l 

read(io)  iz(n,PU),u(n,RH(i),s(i),A(i),isi,in 
I.6yM.l?JR(  U 

I0UM2*1R(2) 

IDU^lxCZ(IL.UMn.zt2))/FL0AT(IUUM1.2) 

Z0UM2r(Z(i0u  42)-Z(  IOUMD) /FLOAT  (Il)UM2-IUUMi*l) 

OZMXNsOZMJiMtAMlNl  (ZPU^liZ0UM21*Qi'»9R 


DZMAXxDZMAXaDZMINai.OOI 
GO  TO  9 

1 FORMAT ( 1RI5) 

2 FORMAT (6E 1 0, 0 ) 

3 FORMAT (215rRE10,0) 

% ^ rnr-,  i r a % 

22__FpRMAT(_RE  1 5^  0 ) 

203  F URMAT ( 8E 1 2, R ) 

1632  FORMA_T(1_M1/////1RX*m  u_jf_Z  L E __SJ.  A ST  ANA  L V_  S I S*//) 

200  FORMAT ( 9Xf *ST  ART 1NG  STEP  *«IR//9X» *F INAL  STEP  x*lR//9X*PKINT 

I INTERVAL  «*IR//9X»STARTING  TIME  «aE  1 3 .5  # 2X # AR , 

2//9X*F INAL  TIME  **E1 3.5»2X» AR,//9X 

3*LA3T  PT,  REGION  NO,  1 x*IR//9X,*LAST  PTt  REGION  NO, _2  s*lU// 

R 9XAM0VING  COORDINATE  SYSTtM  OPTION  x*Ib/J 

201  FORMAT  ( 9X*RAD I uS  OF  TriE  JET  xa£  l 3, 5#  2* , A3f //9x 

l*FLOfc  INDEX  x*I2//9X*OR1GIN  OF  SPHERICAL  SYSTEM  s*EI0.i,*  TIMES  * 

LMtZJ 

202  FORMAT ( 9X » AFHEE  STREAM  PRESSURE  xa 1 1 3. 5# 2X , A9 , ✓ /9x , 

I»FREE_STREAM  VELOCITY  **El3,5f2Xf_A6# //9X,aFREE  STREAM  SPEED  OF 

1 SOUND  s*E13,S,2X,Ab,/) 

205  FURMAT(lHi///iix*MUZZLE  EXIT  C0NDITI0nS*I9X*PLUME  CEnTERLHE  maCh 
1 NUMBER*/) 

20R  FORMAT  (5X»T  I Mg  »t>x  aPRESSuRE  *RX a SOUND  SP»6XaHACHa16Xax/Ra9xaMACHa  ) 

206  FORMAT  (9X*  SPEC  IF  IC  HEAT' RATIO  (REGION  1)'  sa>  R,2, //9xaSPECIF  IC  HEAT 

_ 1 P A T 1 0 _I_R E G 1 0 N 2 )_ ** FR  , 2 , / / , 9 X , * S P E C_IF IC  HEAT  AT  CONSTANT  .PRESSURE 

2 (REG ION  1 ) xa,F7,1,2X,A10»A8, 

3 _ //t9X*SPECIFIC  HEAT_AT  CONSTANT  PRESSURE(REG10N 

R2)  xA|F7(1,2x#A10,A8»/) 

207  FORMAT!  9X.ABAKREL  EXIT  PRESSURE  x*EI 3.5#2X, A9, //9X.*BAR R EL 

1 EXIT  MACH  NO,  x*E13,5//9XaBARREL  EXIT  SPEED  OF  SOUND  x*El3.5,2X 
2# A6)  __  _ 

500  FORMAT(RE12,R«ex,2E12.R) 
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501  F0RMAT(56X,2U2,«) 

502  FORMAT (4E 12, 4) 

. END 

TM  184 


SUBROUTINE  INT  PT 

COMMON/A/ZC200)_#P(200),»i(200)#RM(200b3(20p)#_A(200),»»(3) 

COMMON/B/ZN(206)  #PN(200),UN(20Q),rtHN(2CO),SN(20b),AN(200)#  «N(3) 

COP^ON/C/Z 1(200) , PI (200) ,UI(200),RHI (2QQ) , SI (200) ,AI (200) 

CdM*bN/D/IR(2),GAP,xJ,CP,6Z#0T 

CQPP0N/E/8T  § AL 

COMHON/GA/  GAMl , GAM2, CPI , CP2 

CQMPQN/ZZERQ/ZZ  ERG,  1M0VE#0Z  ZERO»DTZZER,FA  CJ[ 

AL«,5 

BTs , 5 

QAMsGAMl 

CPsCPl 

NCL=2 

NCM=lR(i) 

88_CONTINUt_  

NCll*NCLtl 

NCM1SNCW«1 

MSNCM41 

DO  20  I=MCL »NCM ; 

U I (I ) s U (1) 

P l(l)sP  (I)  __ 

S I(I)=S  (I) 

A I(I)=A  (I) 

20  NHI(1)bKH(1) 

„ DX*1. /FLOAT  (N'CM^iCLJ 

LOOP*0 

DEl*ZN(NO)«ZMNCl) 

DD«1 ,/DEL 

8 XlrQ, 

IF (LOOP  .tU.l)  GO  TO  30 

kNKeMK)  

hNKPtfi  (K+  1 ) 

GO  TO  31 

30  nNKsftMK) 

hNKPswn(K*1  ) 

31  CONTINUE 

DO  1 I*NCL1,NCM1  

XlcxltDX 

NM1*1»LOOP 

NPl*NPltl 

EEeDD*£XP(PI(1))/RHI(I) 

CC*00*(U1(I)*(XI-1,)**NK  -Xl*rtNKP  ) 

FF=0D*GAM  

' PX*(Pl(NPt)-Pi(NMlJ)/OX 
UX*(UI(NPl)-Ul(ivMn)/OX 
SX=(S1(NP1)-S1(NM1))/0X 

W1«Z(1) 

IF (IUOP.GT ,0)  N1«ZN(I) 

RlcRl«ZZERO  

U R E L 1 s’u  1 ( I ) - F 1 0 A T ( I M 0 V E ) * D Z Z E R 0 

PTs-(CC*PXfFF>UxtXj*GAM*URELi/Rl) 1 

UTs-(CC*UX*EE*PX ) 

ST*-CC*Sx 

ZN(I)sZN(NCL)tXI*DEL 

IFUOOP.EG.l)  GO  JO  _7 

PN(l)sPl(I)tP1*DT 


I 

i 

i 

« 
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TH  184 


aSTdfllJII JiiaSr? 

CKSALUG(HH  (I))-*-  UJ/GA M 

B H N L 1 J = P N ( 1 ) / G A h L-  C S N LL)  ? S_i  U ) / C P.»  C K 

RHN(I)=EXP(RMM1)  ) 

AMl)  = SQRJ(GAM*EAPiP'giiij/ftH^I))  

GO  TO  2 

1 -EMm. 5«  tf  1(1  HP  m+PT«DT) 

UN(I)  = .5*(UI(l)fU(muT*0T) 

C SN  m=.5v(  SH  !).♦$(  IJiiMI ) 

10=  1 

I F.(.LLN.(  U_*  U T » 0 , ) 1 0 = - L 

CALL  LPOINT(1,0.,1U,M) 

SNtn=Slt») 

CKCALUG(KH  (I))-P  (IJ/GAM 

RHN(1)s P N c I ) / G A m- ( 5 N (. J J - S J in/ CP ♦ c K_ 

RHN(1)=EXP(RMN(I)) 

ANU^=SOHI(GAM*tXPCPNinj/'WMN(UJ 

2 CONTINUE 

LXQMiMti 

IF(LOOP.EQ.l)  GU  10  10 

DO  l UNCL.NCM  ___  _____ 

U I(I)=U  N ( 1 ) 

F_J.OlSP-_.MU) 

S J(I)=S  N ( I ) 

A .. I ( I ) s A N(I) 

RHl(I)sRHMl) 

J COMUNUE 

L00P= 1 

iXLIJA  e 

10  1F(K.EQ.2)  GO  TU  66 

KSK ♦ 1 

NCLsNCM+1 

NCM=JR(2i  _ 

GAMSGAM2 

CP  = CP2  

GO  TO  88 

6b  RETURN 

END 


TM  1 84 

SUBROUTINE.  LPOINT  ( N , OPT , 10, H) 

MMMflMV  11 20 U 12M1 # HH1 20 Ql/5( iQPlj  A ( g o 0J , ft c 3 1 

COMMON/B/ZN(200),PN(200),UN(200),RHN(200),SN(200),AN(200)»«N(i) 
CflMMjQN/.C/Zl  (iOO^PJ  (20Q)«U1(200),RH 1(200), SI(20U),AI(2QQ) 

CO*MON/U/IR(2),GAM,xJ,CP,DZ#DT 

COMMQN/E/BT*  al 

KJ*| 

1F.CQPU  E <3j_0.  

FIOsFLOAT(IO) 

I*N 

« KSI-IO 

U=HT*(UN(I)»OPT*AN(in 

tMlRsAL*(UCI)*OPT»A( I ) JfEl 

_ EH1L«AL»(U(K)»UPT»A(R))»E1 

ZBs(Z(I)+Z(*))/2. 

*1P*CL _ 

68  CONTINUE 

PAT«CZB-Z(l))/(Z(K)-ZCin 

IF(ABS(RAT).l£,l.)  go  to  1 

>!*ITE(6.  111) 

111  F0RMAT(5X,*CHARACTERISTIC  SHOT  BACK  IS  OUT  OF  RANGE  IN  SUBROUTINE 

1LP0INT*1  

•■■RITE  (6,1  13)  I,K,  IU»ZBf  Z(K)f  Z(I) 

1 13  FORMAT C It  *j  10#  Z6, Z (* ) ,Z(  1 ) *315,,  3E13.A) 

1 EP1*EM1R*RAT*(EM1L-EH1H) 



ZBsZN(I)«EMl*DT 

IFJKJ.GI.OJ  GO  TO ,2 

IF(FIU*ZlI)-FIO*ZBJ  3,2,2 

UUil 

IOs-10 

GO  TO _« 

2 CONTINUE 

It S 3 I 8 A 8 S ( ( Z _B  - 1 B T ) a Z ( 1 1-2LKJ  JJ 

IFITESTZ.LT. .01)  GO  TO  86 

KlPsMP+1 

IF(KIP.LE.IS)  GO  TU  88 

ftR.I  T E ( 6 , 1 1 2 ) 

112  FORMAT (5X, *100  MANY  ITERATIONS  IN  SUBROUTINE  LPOINT*) 

STOP 

86  CONTINUE 

Z1(F)8ZB 

U I(M)8U  C I ) ♦MAT* CU  (K)-U  (13  3 

_ _ _ 3I(M)8S<I)*WAT*(S(K)-S(I3)  

P I(M)sP  ( I ) ♦HAT •(P  (K)-P  (I)) 

CksALQG(RH(1  ) )-P( I )/GAH  

Rh1(M)spj(M)/GA»-(SI(M)-S(I))/CP*Ck 

RHI (M)sEXPCRHI (M) ) 

AICM)sSQRT(GAM*EXP(PICK)  )/RHl(M)) 

RETURN  

END 


TH  I8*t 

SUBROUTINE  MUZZLE C T IM£ , P£ , AE , XME ) 

„.tMMO»4ZJNIT/PM(25).,Avi(25JjXMMt2b)#.TM(25) 

COMMON/FIT/EMP(Z5),EMA(25) ,EMm(25)» INUM 

_j^i(^ioNj:paj#ppijix£Pi?(.iiJfP^ppjtij 

MAX=1 

TD(  nsl  IMF 

CALL  SPLINT (TV, PM, I MUM, TD# MAX , PD, DPD, D2P0, EMP) 

PEsPDClJ 

CALL  SPLINT(TM,AM,INUM, TO , MAX , PD , DPD , 02P0, EMA ) 

._Af*?eii) 

CALL  SPLINT (TM,XMM, INUM,TO#MAX,PO,OPD,DZPP,fcMM) 

XMfcsEfiLLLI 

RETURN 

END 


i 


SUBROUTINE  PUNCH 

COUPON/ A LI  {?  00)  , P (_2  (L01/.U.Ii  0J)J  0)  .S(200)  . *(200)  ,H  L?  1 

COMMON/ 0/ IH ( 2 ) »l»AM#xj»CP»DZ#DT 

CQHPON/GA/  GAMltCAH2jCPI»CP2  

COMMON/ZZERO/ZZERO,  IM0VE,DZZER0.UTZZER#FACT 

I^RITE  (12)  IR(l)»IR(2)»GAMi,GAv*2#CPl«CP2 

HRI TE  ( 12)  (a(N),N=1,3) 

MJL111121  ZliROj  DZZEROfOTZZER 

I1«IR(2)M 

J?R  1 T E C12 )_.UU  ) » P CU#  Ul(J  ) # RHiIJ  iS  ( 1), M U / ill/  1 1J 

RETURN 


END 


SUBROUTINE  SUN 

_ _ SSmQNlt/l  _<  2 o 0 ) j,  P 12.0  OJj  U ( 2.0  0.J  # R H.(  2 0 o ) # S c ? 0 o ) , A ( 2 0 o J J « c i ) 

COMhON/B/ZfJ(2o0)*PN(200)#UN(200),RHN(200),SN(200} » AN(200)#«N(  J) 

__£flKMflN/_p/lR ( 2)  , GAM#  Xjt„C£jJ5Zx01 

Il=IH(2)fl 

PU  1,  1=1,11 

A N(I)=A  (I) 

_P_M.U=P  Cl) 

U N(I)SU  (1) 

S NU)=S  Cl)  _ _ __  

RhN(I)rPH(n 

1 Z N(I)  = Z (I) 

DU  2 1=1,3 

2.  wNa)*HI)_ 

RETURN 
END 


TH  184 

SUBROUTINE  SHOCK ( K ) ^ 

COMMON/  A /Z(  200  )XPJ20  OJ  j UJ200  ) j RM  (200),  S_(  200  ),±(20Q)  fw(  J ) 

COMMON/B/ZN(200),PN(200),UN(200) ,HHN(200),SN(200), AN(200),«N(i) 

COMMQN/C/ZI (200),Pl(200),Ul(200),RHl(200),31(200),AI(200) 

COMMON/ 0/ 1 H (2) » GAM, X J, CP, OZ#  DT 

COMMON/E/BTt Al 

COMMON/GA/  GAM1,GAM2,CP1#CP2 

_C  0 MMON/ZZERO/ZZEROtlMO  VE.U Z ZERO,OTZZER,FA CT 

COMMON /DEBUG/ I DEBUG 

COMMON/ BL AS  I / IBa , B*.C  UN  f SPEED*  R Jfc  T_2_ 

REAL  MHEL#MHEU»MREL2 

MM 

M»J 

IMK.EQ.3)  GO  TO  3 

GAMsQAMl 

CP»CP1 

l«l 

N«2 

OPT«-l.O 

__CO_JO_R 

3 CONTINUE 
GAMxGAM? 

CP*CP2 

1«IR(2)+1 

N*IR(2) 

OPT*1,0 

R CONTINUE 

G2«GAM»I. 

G1«GAM*1 , 

AL« 1 . 

BT*0. 

WN(K)aWJK_) 

MREL* ( U ( I ) -w (K) ) / A ( I ) 

IO*OPT 

LBN- IF  Ik (OPT) 

IPLAGaQ 

STEP=0,r2 

IFLGM*0  

2 CONTINUE 

IF(ABS(MHEL) (LE«1«0vAND,IT(GT(2)  IFLGMal 

IFdFLGM.EQ,  J)  MRELsSIGNU  ,0# MREL) 

ZN(I  )«Z(I)*(ALo/»(N)»BT«V»N(K);>Dr 

1 ZMN)sZN(I) 

_IF(k,EQ.3)  GO  TO  6 

CALL  CL ( I T 1 # BT) 

WN(K)=UN( I )-MHEL*AN(I) 

GO  TO  5 

6 CONTINUE 

CALL  F S 

IF(IBrt.EQ.O)  wN(K)aUN(I)-MREL»AN(I) 

ZDUM«ZN(N) /R JET2 

IF(IBn,EG,l)ir;N(K) >SPEED/SQRT ( 1 • O-BwCON* (SORT ( AlOG( ZDUM  ) )/ZDUM  ) ) 
IF(IbN.EO.l.AND.HT.EO.O.O)  GO  TO  6 
S CONTINUE 

V1*UN(I)-WN(K) 

RH1*RMN(  I ) 

PIBPN(I) 


TM  184 


XMlsVl/ANCI) 

XM12*XMI*XM1  

V12*G2*XM12/(G1*XM12*2.) 

V2«V1/V12 

RH2*RH1*V12 

P2«AL0G((2.*GAM*XM12.61)/G2)*P1 

U2*V2*WN(K) 

AP2=EXP(P2)  

IF(IFLGM.EU.l)  GO  TO  7 

JFl.lM.EQtl,A«.D.h#tG.3)  GO  TO  7 

CALL  LPOINT(N,GPT,IO,M) 

A6«GAM/(AL*AI(1)»HT»AN(N)  ) 

IF(IMOVE.EQ.O)  GO  TO  2l)0 

_AJ*GAM*xj*(AL*(ui(U-027tR0 )/<ZI  (l)-ZZERQ) 

1 /bT*(UN(N)-DZZERO  )/<ZN(N)*ZZERO) ) 

AJ7  = A?*J)J 

GO  TO  210 

2QQ A7«XJ*GAM*QT»(AL*UI(1)/(Z1(1  WZERU)*BT*UN(N)/tZN(  N±rll 

1EHOJ) 

210  _P2_S«P_L(U-Q  p I * A 6 * (U2-UK11)  

P2S=P2S-A7 

AP^S«EXP(P?SJ _ 

ERROR* (AP2-AP2S)/(AP2+AP2S) *2,0 

IF(I0EBUG,EO,2)  wRIT£(6,10)  1T,I,N 

IF  (IDE8UG.E&I.2)  WRITE  (6, 20)  MREL,  ZN( I ) , ZN ( N)  »«N (K)  , UN(  1 ) * AN  ( I ) , 

1V1»RH1,P1,XM1 _ 

IF ( IDEBUG.EQ , 2)  WRITE (6, 20)  XM 1 2 , V 1 2, V2, RM2 , P2, V2, AP2,  Ab , A 7 , P2S 

JFI  [DEBUG, EU, 2)  wRITEIA,20)  AP2S, ERROR, Xj 

IF(AbS(ERROR),LT,t.E-Ofl)  GO  TO  7 

S1*SIGN(1., ERROR) 

S$TEP*SIGN(  1 .,STEP) 

IFtABS(MREL) .LE. 1 . 05. AND.  ABS( STEP ) . NE . 0 . 0 1 ) STEPsO, 0025*SSTEP 
I T*1 T* 1 

IFC  IT .GT .75)  GO  TU_U0 

IF(IT,GT,2)  GO  TO  40 

50  CONTINUE 

MRELl=MREL 

ER 1 sERRQR  _ 

S2*S1 

_ _MREL*MREC_*STEP _ 

GO  TO  2 

«0  CONTINUE 

IFCSl.ME.S2.UR.IFLAG,EQ.l)  GO  TO  45 

IF  (ABS(ERROR),LE, ABS(ERl))  GO  TO  50 

ST EP**5TEP 
GO  TO  50 

110  CONTINUE 

WR1TE(6,111)  X,WREL,U2 ,RH2, P2 

111  FORMAT (5X»  *TUU  MAM  ITERATIONS  IN  SUBROUTINE  SHUCK* , 2X , *K  = * , i 5, 2X , 

J * MREL**, El 3 l5 ,2X,»U  2 **,  E \_i , 5,2  X L*  RH2*  * ,£ 13 , 5, 2 X , * P 2 s * , t 13 , 5 ) 

STOP 

120  CONTINUE  _ 

WRITE (4,121 ) 

121  FORMAT C5X#  * ABSOLUTE  VALUE  OF  MACH  NUMBER  IS  LESS  THAN  ONE*) 

CALL  EXIT 

45_CONT  JNUE 

MREL2*MHEL 1 »ER1 * (MREL*MREL1 j / (ERR0R»ER1 ) 
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1 


TM  18* 

MRtll*MR£L 

ERUERROR 

MR£L*MREL2 

IFLAG«1 

GO  TO  2 

7 UN(N)«U2 

RMN(N)*RH2 

PN(N)sP2  

8N(N)sSN(  I)fCR*(  (PN(N)-PN(I))/GAM.ALOG(RMMN)/HHNCI))) 

_ AN(N)sSQRT(GAM*AP2/i?H2) 

IMBT.GT.O.)  GO  IU  9 

6 CONTINUE 

BT«,5 

ALr, 5 

I T*  1 

1FLAG * 0 

STfcP«0,02 

lFLGMsp 

GO  TO  2 

9 RETURN 

10  FORMAT (1415) 

20  FORMAT  ( 1QE 1 3,5) 

' END 


SUBROUTINE  SHLINt(X, Y#N,EM) 

_OI_ME_NSION  X(N),Y(N) 

DIMENSION  SH(25),G(25)#EM(25) 

DAT*  SB ( 1 ) » G C 1 )/-0.5#0.0/ 

NO=N-l  " 

If (N0.LT.2)  uO  TO  20 
DO  10  1=2# NO 

A « (»OHii»m/ht  

c = (x(i*n»xm  j/6. 

* = 2,'*(AtC)»A*SH(i-l) 

SB( I ) = C/ft 

F = ( Y ( I ♦ 1 ) - Y (l)?/(X(I»l)-X(I))»(Y(I)-Y(I-l))/tX(l)«X(I»l)) 
GU)  s (F-A*U(  1-1 ) )/a 
bM(N)  * G(N-1 )/(2,+SB(N-l)) 

DO  30  I = 2 # N 

K S N+l-I 

tM(K)  = G(K)«SdCK)*£vfcKii)  “ 

RETURN 

END 


SUBROUTINE  SPLINT (X,Y,N,Z,MAX,YInT,DY0X,DY2DX,Em) 

_OIMENSION_  X(N), Y(N),Z(MA  X ) , Y1NT (MAX  )'OYOX(MA  X)#D  Y2DX  (MA  X ) 

DIMENSION  EM(Zb) 

DATA  SRw/0/ 

INTEGER  SRw 

III  « SR* 

DO  1«0  1*1 »MAX 
K = 

IF (Z ( 1 )-X ( 1 ) ) 70 ' 60 ' 90 

YjNT(I)sY(l)  

SK  = X(K)-X(K-1) 

GO  TO  130 

IF(Z(I),GE.(1.1*X(1)-.1*X(2)))  GO  TO  120 

WR I T E_ ( 6,1 OOO)Z(I) 

SRY.  s 16 
GO  TO  120 

KS 

1F(Z(I).L£,(1«1»X(N)-.1«X(N~1)))  GO  TO  120 
WRITE  (6,1000)  2(1) 

SRw  * 16 
GO 
If 
YI 
SK 

GO  TO  130 

K*K  ♦ 1 _ 

IF (K-N)  90,90, BO 

-CONTINUE _ 

SK  « X(K)-X(K-1) 

Y I NT  C 1 3 * tM(K-l)«(X(*  )-Z(I))*»3/6./SK  »EM(K)»(Z(D»X(<-1  ) )»*3/6 

♦/SK+(Y(K)/SK  -EM(K)*SK  /6,)*(Z(I)»X(K-1))^(Y(K-1)/SK  -E*(K»1) 

♦ *SK/6,)*(X(K)-Z(I))  .. 

DYDX(I)*-EM(K-1  )*(X(K)-Zd))»*2/2.0/SK  «ttt(K) * ( X(K*l)-Z ( I ) ) **2/1. 

♦ /SJ'  ♦ < vj  K_3  -LV  c K • l > J / SK  -(EM(K)-EM(K-l))*SK/6, 

" D Y2 DX ( I)ViM( K*1 ) * (X ( k ) - Z (I ) ) / SK ♦ E M ( K ) * ( Z C I ) - X ( K - 1 ) ) / S < 

CONTINUE 

MX A s MAX 0 ( N , M Ax ) 

IF(SMw,EU,16)  WRITE (6,  1010)  N , M A X X ( I ) , Y H )_»  Z ( I ) , Y I N T ( I ) , D Y D X ( I ) , 
♦DY2DX (I ) , 1*1 , MXA) 

SRw  * III 

RETURN 


1000  FORMAT  (59H  SPLINT  USED  FUR  EXTRAPOLATION,  EXTRAPOLATED  VALUE  = , 
♦614,6) 

1010  FORMAT  (2X,21HN0,  OF  POINTS  GIVEN  «'I3'iOH,  NO,  OF  INTERPOLATED  PO 
♦ INTS  ’*,13/1  OX*,  1MX,  19X,  JHY,  16X,  1 1 HX  - 1 N T E RPOL , , 9 X , 1 1 HY-I  mT  fc  WPOL  . , 

♦SX , 14MD VOX»I NTtRPOL , , 6X, l5HDY20X»INTfeRP0L . / ( 6E20 ,B) ) 

END 


TH  184 


SUBROUTINE  STEP 

jcj^mun/a/z  (20  o)j  p c2  o qj  ^ u (200 i jL.m  2mi,  su^AjAHslouaL}.! 

C0PM0N/D/IR(2) #GAM, XJ,CP,0Z,3T 

CQMPQNZ6/KK#.U.#K0,  TI*E*Jl*fc£ 

COMMON/DZ/  DZ^I^UZMAX 

CQMMQN/GA/  GAMi ,bAMg.CP1  .CP2 

COPMOM/l)TSTAb/OTSTA3 

DL=L.ef06. 

KS1 

NCL=2  _ _ _ 

MCMs 1 H ( 1 ) 

5 CONTINUE 

DZs(Z(NCM)-Z(NCU  )/PLOAT(NCM.NCL) 

IF  (DZ.GT.DZM1N)  GH_TO  10 

IF(NCL.NE,2)  GO  TO  20 
IF(IR(1J.LE.«)  GO  T.0  JO  

GAKsQAM J 

C£=£P1 

CALL  INTER(i#lR(l),-l) 

mu  = iR(u-j 

NCP=IK(l) 

GO  TO  5 

20  GAMsGAF'2 

C.P.g.LP.2 

IF(IR(2)-IR(1).LE.3)  GO  TU  10 

....  CAUL__IKIEP(IK(U,IRC2Jj^U  ... 

NCM*IR(2) 

_eo  _iu_  5 

10  CONTINUE 

1 F1P Z « LlJ?Xl*x j_. (,0- _LQ.  go 

IF(IH(2).GE.199)  GO  TO  «0 

JMNCL.NE.2),  GO  TU  10 

GAMsGAMl 

CPsCPJ 

CALL  INTERUf  lKCD'M) 

IW(1)«IH(1)»1 

NCMsIHd  ) 

GO  TO  5 _ 

30  G AMsG AM2 
CP*CP2 

call  INTER ( 1 R ( 1 ) t IH(2)#*l) 

N<>=  I R ( j?J 

GO  TO  5 

«0  CONTINUE  __  _ 

DO  7 IsNCL  #NCM 
Ul=  _ U ( I ) _ 

A 1 2 A ( 1 ) 

XNlt A H S(Ul / A 1 ) ♦ 1 . 

DT1=0Z/A1/X^1 

_ IFiDTj.LT, DTJ  J>T*DH 

7 CONTINUE 

IF(K,EQ,2)  GU  TO  6 
If  (IH(2)-IR(1),LE,2)  GO  TO  8 

K«K»1 

NCL=NCM+ J 

_NCMsIR(2I 

GO  TO  5 


RETURN 


IM-lfllt 

SUBROUT  INfc  IM  ER  ( 1 1 , 1 F , IOPT  ) 

MMF'QN/A/2(200),PUOO)JUt200]#(:!H4?_O.Q)/S(200)/A(290Jj!N(3) 

COMMON/ d/ZN( 200) , PM  200) , UN (200 ) , RHN(  2 00 ) , SN(2Q0)  , AN(200),*N(  3) 

C£)mmj)N/d/IR<2)  »GAm,  xJ,CR,OZ,OT 

COMMON/F/RJ£T,XiVi£,RE,  AE,PINF,  A I *F , U INF , RH1NF 

Illsli  + l - 

1 12=  1 1 12 

1FN=IF+I0PJ _ _ 

J I = 1 1 1 

l DUM=  1 

IF(IF.tU.IR(2))  lCU*t=IR(2)tl 

DO  1 1«112,IFN 

RATsFLUAT  (1-H1) /FLOAT  (1FN-II 1 ) 

ZNmiZ (J J 1 ) ♦ ( z u F ) ri  LLL1JL1_?MT 

DO  10  JsJI,  IF 

_ JusJ__  _ 

IF  (Z(J)-ZMl  ))  10,3,2 

10  COM1NUE 

IF  (Z(IF ) ♦ , 000 1 , G 1 • ZN  C 1 ) ) GO  TO  3 

wRI  Tt  16#  1 00  ) _ _ 

100  FORMA T ( * ERROR  IN  SUBROUTINE  INTER* ) 

STOP ___ 

3 U N(l)su  (JU) 

S N(l)sS  (JU) 

P N(I)=P  (JU) 

A N(J)  = A (JU)  _____ 

RHN(I)sRH( JU) 

JlcJU 

GO  TO  1 

— 

HATs(ZMI  W(Jl))/(Z(  JU)-Z(JD) 

(INXIJSU(JL)*(U(JU)-U(JL))*RAT  __  _ _ 

RHN(I)sRH(JL)*(Rm(  JJ)-RM(JU)*HAT 

__  P2*tXP(P( JL)  ) ____ 

PN(i)=P2f(EXP(P(Jij))-P2)*RAT 

SN(1)«S(IUUM)*CP*(  (ALOG(PNU)  )-P(  IDUM)  ) /G AM-ALOG ( RHN ( I ) /RH ( I DUM ) ) ) 
AN(I)sSQRT  (GAM*Pn( I ) /RHN ( I ) ) 

PN(I)=ALOG(PN(I)) 

J1=JU 

1 COM  INUE__ 

,MslFM 

IR21slR( 2)  ♦ 1 

DO  U 1 = 1F 1 , 1R2  1 

J*I*IOPT  _ _ 

Z N(J)=Z  (1) 

U MJ)=U  (I) 

S N(J)*S  (I) 

P MJ)sP  (1  ) 

A N ( J ) = A (1) 

U RhN ( J ) sRM ( I ) 

I R ( 2 j * I W ( 2 ) ♦ 1 UP  T 

IR21=1R(2)*1  _ _ ' 

DO  5 1*1 12, I R2 1 

Z (l)sZ  NT  1) 

U (I)su  N ( I ) 

SJJi*S  N ( I ) _ ; 

P (I)*P  N ( I ) 


U?  = XM*»«A2 

DttT  = TINIT-JM(l) 

iterz=i 

ZLOV.sZE 

ZUPsZtf  U1*T1NII 

CONTINUE 

ZMDsZLOr.*0.2*(ZUH-ZLOw) 


DETERMINE  MACH  NUMbER  AT  MACH  DISC-POluT  1 
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i 


I 

TH  1 8*1 



ZM(  1 )=Z*D/HJt  T ♦ 1 , o 

CALL  SPLINT  (ZZ#Z<^CL,  INUMP,  ZM,  MAX  , XM,  DXM,  D2XM,  fc^CL) 

X_MMD=XM(  i J 

Cr*MMM*«M*0*«*tM*00**0**«M«**t**MMM*t»MM*0**»MMM**«M 
C LOCATE  ORIGIN  OF  STREAMLINE  A?  MUZZLE  £XIT 

ITEPSsl 
RAT  10=0 ,5 
20  CONTINUE 

TDUPS=RAT1U*TINIT*U  . 0-RAT  10)*T1 
UE=RATIU*U2+( 1.0-HATIO)*Ui 
TUUNP=T  IMT-IZMO-ZD/UE 
ERRsCTDU^S-TDUMPJ/TDUMS 

IFCTUUMP.lT.T1, Ok. TDUMP.GT. UNIT)  GO  TO  650 
IF(ABS(EHR),LT.fcRHVAL)  GO  TO  300 
ITfcRSsITERSfl 

IFCITEKS.GT.20)  GO  TO  200 

RATIUS(T1-T0UMP)/CT1-TIN1T) 

GO  TO  20 
200  CONTINUE 

IFCIOEBUG.EQ.i)  *RI TE C6, 2000)  

2000  FORMAT ( l X » * T 00  MANY  ITERATIONS  FOR  LOCATION  OF  STREAMLINE  IN  SUBRO 

1UTINE  INIT«) 

CALL  EXIT 
300  CONTINUE 

PE=RATI0*P2M1.-RATI0)*>1 

_ jAE*RA T IO*A2* ( l,-RATIO)*Ai  

XME*RATIO*XM2+(l ,.RAT I0)*XM1 

C DETERMINE  PROPERTIES  AT  POINT  l 

C*** ***********  ********************************************************* 

Fl*l  ,0*0,5*G1*xmE**2 
F2®1 ,0*0.5*GI*xmmd**2 
PPEs(F1/F2)**G3 
RHREs(F1/F2)**Gm 
PCI) sPPE*PE 
RH£sGAMl*PE/AE/AE 
RH( l )sRMRE*RHfc 

A(n*SORT(GAMl*P(l)/RH(I)) 

U(l)«XMMD*All) 

wITERsU(i)«I.IO*A(l) 

I'FCIMOVE.FQ.O)  GO  TO  15 
ZZERCtsZMD*FACT 
DZZEROs NITER* FACT 
DTZZEKsO.O 
15  CONTINUE 

IF(IDEBUG,EQ,3) 

r«RrtE(6,l)  Pt,AE,XME#P(n»»H(l),A(l)#U(l)#**ITO 

c** **** *************** **************************«*********************** 

C 'DETERMINE  PROPERTIES  AT  POINT  2 

C************ **********************************  *********** *«***••• ****** 

XF  PT 1*1,10 

Xls(2,*GAMl»xMPT 1**2-G1 )/G5  

P(2)sxi*P(l) 

FXl=(G5*Xl*Gl)/(Gl*XIfG5)  

U(2)swlTER*(U(l)-wITER)/FXl 
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RH(2)sFXI*Rh(1) 

CONSMO=RH(2)*u(2)*(ZMD-ZZERO)**JM 

IF ( IDEBUG.EU, 3) 

1WR]T£{6,1)  XMPT1,x1,PC2)jFxI,U(2)^Rm(2)  _ _ 

C’ 

c 

DETERMINE  BLAST  *AV£  PROPERTIES 

* 

r , 

XlBwsP(2)/PINF 

FXllr  tG7*XIBrt*ub)/(2.  *GAM2) 

VSBw  = S(JRT(FXp  *GAM2*P1NF/HHINF) 

FXI2*(G7*  XlBn+G6)/(G6*xlB*tG7) 
HH(IR2)=FXI2*RhINF 

UDUM2sVSBk-VSBa/FXI2 

ZSHK*ZJj  + VSBfl*DELT  _ _ 

- 

CONSB*=RH(IR2)*UCUM2*(ZSHK.ZZERO)**JM 
C DETERMINE  PROPERTIES  AT  MO  SIDE  OF  CONTACT 

CM*00M«M0«  *««>  * * ^**  * *_*  ** 

CUNST=C0NSMD*DELT/kM(2) 


__cculq_uo2#  101 ,100)  ,joyM 


100 

_ yv * V.  A A V.fcl  * V » f • v v / IJrf  • 

CONTINUE 

ZC*0,5*(ZM0fZS«K) 

• 

30 

ITEHZCsl 

CONTINUE 

YDUMsZC**3-2,*ZZEHO*ZC**2*ZZERO**2*ZC-CONST 

DYDUMs3.0*ZC**2-<4fO*ZZtRO*ZC*ZZERO**2 

ZCDUMs2C-(>DUM/UYDUM) 

ERRZC*(ZC-ZCD'jM)/2C 

IF { I DEBUG , E0 , A ) WRITE(6,3)  I T ERZC , ZC , ZCOUM, 
IF(ITfRZC.GT.20)  GO  TO  «0 

YDUM,DYDUM, 

ERRZC 

I TERZC* I TfcRZCf 1 

IF(ABSJERRZC) .LE.ERRVAU)  go  to  50 

ZC=ZCDUM 
GO  TO  30 

«o 

CONTINUE 

IF ( IDEBUG.E0.3)  *R I TE ( 6 . 5000  ) 

5000 

FORMA  T( 1 X » *T00  MANY  ITERATIONS  FOR  CONTACT 

unit*) 

POSITION  IN 

SUBROUTINE 

101 

GO  TO  650 
CONTINUE 

ZC=(ZZERO*SURT(ZZtRO**2*«.*CONST) )/2,0 
GO  TO  55 • 

102 

CONTINUE 

ZCsCONST 

50 

CONTINUE 
UC1 =ZC/UEL  T 
PC  1 sP (2) 
RMClsRH(2) 

c 

c 

DETERMINE  PROPERTIES  AT  Bw  SIDE  OF  CONTACT 

************ 

c 

*•**1 

PCZsPCl 

- 

UC2SUC1 

RHC2sC0NS8a/(UC2*IZC-ZZER0)**JM) 

c 

c 

c 

***«i 

CHECK  SOLUTION  USING. DENSI TY  AT  brf  SIDE  OF 

CONTACT 

* * * * t 
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I 


ITERZslTERZM ,n  ,u" 

CRRQR«(RH( I R2 ) »RHCZ)/RH(1R2) 

IF(IDEBUG,EQ,3) 

_ 1KRITE(6»3)  IURZ, * 1 8 w ,_y  SB w » RM ( 1 R 2)  »RHC2fERKQR • ZMQ» ZC » Z Si M < 

IF (A6S( ERROR ) .LE.tRRVAL)  GO  TO  700 

IF(ITERZ.GT.20)  GO  TO  600 

IF (ERROR, GT  ,0,0)  GO  TO  500 

_ ZLOwsZF'D 

GO  TO  10 

500  CONTINUE  

ZUMsZWD 

GO  TQ  IQ 

600  CONTINUE 

_JF(I0EBUG.E3,3)_wRl_TE<6,3000> 

3000  FORMAT (1X,*T0U  Many  ITERATIONS  FOR  SHOCK  VELOCITY  IN  SUBROUTINE  IN 

__  1 1 T») 

650  CONTINUE 

ITERTsITENTM 

IF(ITERT.GT.20J  GO  TO  675 

T I NTT  * I JUMP  *T INI T 

GO  TO  5 

675  CONTINUE 

»RITE(6,40~00) 

4000  F QRWAT ( 1 > » »TUU  MANY  ITERATIONS  IN  SUBROUTINE  INIT  FOR  THE  INITIAL 
ITIMt  STEP* ) 

CALL  EXIT 

700  CONTINUE 

IF(ZC,LT,ZMl).ORfZC.GT,ZSHK_)  GO  TO  600 

IF(ZMO.GE.OMIN)  GO  TO  750 

TlNlTsTJUMP«TINIT 

GO  TO  5 

750  CONTINUE, _ 


C DETERMINE  ALL  FLOW  PROPERTIES  AT  THE  DISCONTINUITIES 

TIMErTINIT 

W( 1 )*WITER 

Z(l)»ZMO_  

S(1)=0.0 

A(l)sSQRT(GAMl*P(i)/RH(l^))_ 

P(l)*ALOG(P(l))  " 

Z(_2)_=ZU_i 

S (2 ) *S ( 1 ) *CP 1 * C ( ALOG ( P (2 ) ) »P ( l ) ) /GAMl-ALUG CHH(2 J /RH ( 1 ) ) ) 

A(2)xSURT(GAM1*P(2)/rh(2^ 

P 1 2 ) * A L 06 ( P ( 2 ) ) 

Z(IR1)*ZC 

S(IR1J*S(2) 

P(IRI)SAL OG(PCl) 

RH ( 1 R 1 ) sMhC 1 

U (I R 1) « U C 1 

A(  IRlJsSORf  (GAM1*PC1/RHCD 

rt  ( 2 ) sUC  1 

Z(IR2P)=ZSHK 

P(IR2P)sAL0G(P1NF) 

U(IR2P)=U1NF 
RH( IR2PJSRH1NF 
S(lR2P)s0.0 
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l 


A(IR2P)rSQRT(GAM2*PINF/RHINF) 

M3)«VSB* 

Z(IR2)s2(lR2P) 

.■■■EL(IP2)«PL£j. 

U(IR2)*U0UM2 

S(lR2)=CP?»((PClR2)-P(IR2Pn/GAM2.flLQGfRH(IR2)/RHlNf))*S(IR2P) 
A(lR2)sSQRTlGAM2*tXP(P(IR2) )/RH(IR2)) 

P(IR1P)=AL0G(PC2) 

__WJRiPl*UC2 

RH(1R1P)*RH(IP2J 


S(IR1P)=S(1W2J • 

A(IRlP)rSG«T(GAM2*PC2/RHC2) 

IPjjMIfJR  ( U : 

IPUM2dR(2) 

___ZPy_MX*(.IUPPyt)-Z(2i)yfi^0fI(I0UMl-2)  ■ 

Z0UV2s(Z(IDU^2)-Z(10UMl ) ) /FLOAT  ( ll>U*2-I  OU*l  • i ) } 

PZHlNaDZMlM»AMlM(2DUMli2DUM2  3«Q.99Q 

DZMAXsOZ*AX*DZMIN*l ,001 

RETURN 

1 FORMATClOt 1 3,b) 

_2_F0RMAU1 9 15 1 

3 FOBMATCi I5,8tli,5) 

END 


c** 

SUBROUTINE  ACOUS(ZCON,Z8rt,ufrdj?Xo,HHORCT) 

c 

c 

VAR(l)  IS  TIME 
VAR(2J  IS  Z 

c** 

common/f/rjet,xme,pe, ae.pinf,ainf,uinf,rminf 

COMMON/G/KK,LL,KO,TlME,TFINAL 
COMMON/G A/GAM l , GAfr>2,  CP  1 * CP2 

\ 

COMMON/BL  AST/ IB*v,8*vC3N, SPEED,  RJET2 
COMMON/KUTTA/  CUVAR(2)f V AR ( 2 ) , C I f DER C2 ) , I I 

- 

DIMENSION  ERRVAL(1),ELT(1),ELE1(1),ELE2(1) 

DIMENSION  RFIT(20),PFIT(20),RHOFIT(20),RFITU(20),UFIT(20J 

DIMENSION  EMFITP(20) ,EMF1TR(20),EVFITU(20) 
DIMENSION  XDUM(1).YDUM(1),0YDUM(1),D2VDUM(1) 

• 

DATA  IF  I T , IDUM/20, 1 / 

DATA  (RFITI I),  1*1, 20J/-50.0, -45. 0f -40.0, -35.0,-30.0,-27. 5, -25.0, 

i 

1 -22. 5, -20. 0, -1 7. 5, -IS. 0, -12. 5, -10. 0,-7, 5, 

2 -5. 0,-2. 5,0. 0,2. 5, 5. 0,7. 5/ 

DATA  (PFIT (I), 1*1, 20 J/0.0, -0,16,-0,24, -0.«2,-0, 62,-0, 78, -0.90, 

1 -1.03, -1.20, -1.53, -1.90, -2. 35, -2. 83, -3. 22, 

2 >3.40, -2,90, -0,60,4.30, 13, 0,23,1/ 

DATA  (RFITU(I),I*l,201/-50.0,-42.5,-3S.O,-27.5,-22.5,-17.5,-15.0, 

1 •12.5,-10(0,-7,S,-5,0,-4,0,-2,5,-l .0,0.50, 

2 2.0, 3, 5, 4,5, 6.0, 7,5/ 

DATA  (UFIT(l), 1*1,201/0.0,0.04,0.095, 0.162, 0.22, 0,30, 0.35, 0.41, 
1 0.52,0.65,0.77,0.825,0.91,0.97,1.0,0.955, 

2 0.755,0.55,0.045,-0.725/ 

DATA  H0L1/2H  / , HQL2/2HB*/ , MOL3/2HMD/, H0L4/2H  C/ 

EXTERNAL  DERSUB, CHSJB 
DO  8 I* 1 , I F I T 

RH0FIT(I)*PF1TII)/GAM2 
8 CONTINUE 

CALL  SPLINE(RFIT,PFIT,IFIT,EMFITP) 
CALL  SPLINE(RFIT,RmOFIT, IFIT,EMFITR) 

CALL  SPLINE (RF ITU, UF IT,IFIT,EMFITU) 
IMD*0 

N*1 

ITEXTsO 

NT*1 

ELT ( 1 )*TF INAL 

SPEC *0.0 
C I*DT 

CIMAXsO.O 
ELEl ( 1 )*0 • 00 1 

ELE2(1)*0.0 
VAR(i)*TIME 
VAR (2 j*ZBh 
11*0 

- — — — 

3 CONTINUE 

CALL  INT1A(H,N,NT,CI,SPEC,CIMAX,  I ERR,  V AR, CUVAH , PER,  ELE 1 , ELE2 , EL T , 
1ERHVAL, DERSUB, ChSUB,  I TEXT ) 

KO*KO«l 

GO  TU( 1 , 1 , 2, 2) , 1 ERR 
1 CONTINUE 

IHK0,GE,*K,0R,VAR(1).GE.TFINAL>  &0  TO  5 
1F((K0/LL)*LL.NE,K0)  GO  TO  150 

5 CONTINUE 
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*RITE<6# 1000)  KO»VAR(l) 


Z0UMsVAR(2) 


■ IIII.li 


CALCULATION  OF  T Ht  MACH  DISC  LOCATION 


lF(IMD.EO.l)  00  TU  b 

CALL_MUZZLE  ( TOUVj  PEX  AL*  XMJEJ 

PRAT  sPi/P I NF 

LF(PRATjGE,0^0)_G0_T0_7 

IMD*  1 
GO  TO 


7 CONTINUE 

ZPO*CONMD*XVfc*SOHT (GAMJaPHAT) _ 

IF(ZMD.LE.UO)  I.*'Osl 
6 CONTINUE 


WR I It  16. « Q.QAL2 

IF(IMD.EQ.O)  wRI Tt (6# 1 0 ) H0L3#ZMD 

_ HOLpHULa . 

00  100  1*1,20 


!(<Q]X'ElIlulr4UI<] 


RDUM=ZCON  + RATIO*(ZOU'1-ZCON) 

_XOUP>U)*iROUM-SPEEO*TOUMl/(R,IET?<LALOGiRDU'</RJtT2i) 
IF  ( XDUM(  l).GE.RFIT(l))  GO  TO  2*0 

PRESS*!  iO 

RMOsl.O 


GO  TO  30 

20  CONTINUE  _____  

CALL  SPLINT(RFIT,PFIT,IFIT#XDUM,  IDUM,  YOU*# D YDUM# D2 YDUM# EMF I TP ) 

_ PRESSsYOUM( 1 )*RJET2/R0UMtJ ,0  __  

CALL  SPLINT (RF IT »RHOF IT# IFIT»XDUM, IDUM#  YDUM# DYDUM, 02 YDUM, EMF I TR ) 

RHQ»Y0UM(l)*RJET2/RDUM»l  .0  

CALL  SPLlNT(RHTU#UFIT,IFIT#XOUM,IOUM#  Y0UM#DY0UM#02Y0U^,EMFI  TU) 

UVAL*RM0-1  t0-YDU*(U*lRJET2/RDUMJ**2*RHURCT_ 

30  CONTINUE 

IF(I,EO,20)  H0LSH0L2 

ZDUMMY*HDUM/RJET2 

wRITE(b#10)  HUL,Z0UMV1Y,  PRESS#  RHQ.UVAL . 


H0L*HUL1 

100  CONTINUE  ; 

if(ko.ge.kk,or.var(i),ge.tfinal)  go  to  a 

ISO  CONTINUE  

GO  TO(3#3#4#«)#II 
C0NT1N 


WRITE (6»2000)  IERR 

CALL  EXIT  

a CONTINUE 

wRITt(6_*3000) \ 

STOP 

_ . 10  F0RMAT(lX,A2,aE13.S) 

1000  FORMATdHl#  10X#*STEP  s*  # I « # 1 0 X , * T I ME  **#E13.S#//) 

2000  FORMAT (//# 1 x# 26KI NT I_NON •CONVERGENCE# I ERR* #12) 

3000  FORMAT!//# IX, 17HtNO  CONDITION  Mt T ) 
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TM  18A 

SUBROUTINES  1A(  1 I , N, NT ,C I , SPEC , C IMAX , IERR  , VAR,CUVAR,DER,  ELE  1 , 

lELE2,ELT,ERRVAL,DEPSUB,CHSUB,  I TEXT)  _ 

*********  DOCUMENT  DATE  09-01-68  SUBROUTINE  REVISED  08-01-68  •**•****< 

DIMENSION  SlVAR(20),SEL£lJZ0),SLEt(20),ELE2(j>0),DER(21), 

1 FDERV(21),S0Y(20),SDY1 (20), V1NCH(20),E*RVAL(20),£RV0VH(20), 

2£LT(10),S£LT(13),RELMIN(20),STEP(3) 

DIMENSION  VAR(21)#CUVAR(21) 

INTEGER  TEX ( lb) 

INTEGER  CODE,  TPSri,  SURHAF , STEP,  TEST, DCODE 

REAL  Kl_ 

C BEGIN  INITIALIZATION 

IF  (II  ,GT,  0)  GO  TQ  520 

1 TP*0 

SSPEC  = SIGN(SPEC,CI ) 

SCIMAXeSI GN( ClMAX , C I ) 

__V_ARJ  s VAR(l)  _____  

IF  (Cl  .EG.  0.0)  GO  TO  151 

IF  (SSPEC  .EG,  0.0)  GU  TO  7 

IF  ( ABS  ( SC  1 MAX  ) , G T , ABS(SSPLC)  .OR.  SCIMAX  ,EU,  0.0~) 

lSCIMAXxSSPEC  

C TEST  TO  SEE  IF  VAR  IS  2ERO 

IF(ABS(VAR1).GT.1.0E-H)GQ  TO  2 

TPsSSPEC 

GO  TO  7 

2 IF  ((V AR1/SSPEC)  , GT , l.E-13)  GO  TO  <4 

_3__K_1*0. 0 _____ 

GOT  06 

«_K1*L,0 _ . 

6 TP*VAR1-AM0G(VARI , SSPEC) 

IF(ABS(TP-VAR1).LT.1.E-12)K1e1.0 

TPeTP  + M *SSPtC 

_I_F  (ABS((TP-VAH1)/VAR1)  ,LT , l.E-11)  TP  * TP  ♦ SPEC 

C TEST  FOR  DIRECTION  OF  INTEGRATION 

7_M*lf0 

IF  (Cl  .IT,  0.0)  KIe-1.0 

CIkeC1*K1 

CIMAxksSCIMAx**! 

TPKtTP*Kl  

VARKEVARl*kT 

C SET  UP  STORAGE  FUR  INTERNAL  USE  

NP1  ENVi 

NELT  = 1 

REMAlNsO.O 

NhAF  *0  _ 

Ni'SiNt  ~ . 

SUMMAFsO 

LGOP*0 

DO  91  1=1,3 

91  STEP(1)*0 

IERR»1 _ 

DO  8 Ui.NPl 

8 CUVAR(I)  svAR(I)  

DO  101  I * l , N 

101  3ELE 1 ( I ) sELE 1 ( I ) 

IF  (NT  ,EO,  0)  GO  TO  13 

100  IF  (NT  #EO,_i_)  GO  TO  10 

NTM1*NT-1 
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Th  m 


EITK*K1*ELT(1) 

PS)  9_1=L«  NTM1 

ELTK2sKl*ELT(Itl) 

If  CELTK  .LT.  ELTK2)  GO  TQ9 

GO  TO  500 

? ,ELTK«EU*2 

10  CONI INUE 

E-LIK.s|<l»fcLT  (NELT) 

If  (VARk  ,LT.  ELTK)  GO  TOll 

if IN E LI., E ftt.wn  GQ  TO  13 

NELT*NELT*1 

GO  TO  10 

11  NEITI  = NT-NEITM 

GQ-1QL2 

13  NELTL*0 

12_iiQ_l<l_I*ltN__  

14  RELMIN(I )=$ELfcl(I)/i28,S 

If  (NT  .EG.  0)  GO  TO  996 

DO  995  Is  1 # NT 

i.95_?imn*ELui ) 

996  CALL  DERSUb  ( 1 1 , CUV AR, DER, N, V AR) 

If_ai„. E Q .„« ) GO  TQ  120 

DO  15  I = 1 , N 

15  FPERV ( I ) spfcR ( I ♦! ) 

11*1 

IE §1*0 

DO  300  I s 1 » 1 5 

jojlje  xj.  j l*  o 

TEX ( 1 )*1 

TEX(2)*i 

KK3*1 

J.F_1I  TEXT!  to  35,.to.3,J>35 

151  wRJTE  (6,1000) 

.IPQlO.  F ohm at_1  /./1LH  _c  i_is_zmi_  __  

STOP 

_c END  Of  INITIALIZATION 

520  11*1 

... TPSH=0 

LTSHsO 

VARK*VAR(1J*K J 

CIKsCUKl 

S1«VARK»CIK 

If  (SSPEC  ,EQ,  0,0)  GO  TO  525 

KKs  1 

If  (NELTL.EQ.  0)  GO  TO  17  * 

IF  (ELTK-TPK)  16,16,17  __  _ 

16  CV*ELT»v 

COOES  l 

GO  TO  16 

17  CV«TPK 

C00E*2 

l*_|f (ABS(CV).LT , 1,E-12)GO^TO_5JO  _ 1 

If  (CV-Sl)20,20, 19 

19  IF(ABS((CV»S1)/CV).GE..1E»11)0Q__T_Q  535 

20  If  (NELTL  .EO.  0)  GO  10  540 
lf(AHS(CELTK-TPK)/Cv),LT,,lE-U)GO  TO  550 
IF  (CODE  ,EU,  1)  GO  TO  545 
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DXSTP-VAR(I) 

JLEX(b):  1 

TP*TP+SSPEC 
TPKsTP*K 1 


TPSHsl 
GO  TO  560 


SHORT  INTERVAL  DUE  TO  BOTH 

_T_P=TP*SSPEC 

TEX (6) s 1 

TPK=TP*K 1 _ 

TPSHsl 

GO  TO  545 

IF  HERE  CV  IS  LIKELY  ZERO 
IF  (SI  ,LT.-l,Ut-12)GO  TO  535 
IF  (CODE  ,EO.  1)  GU  TO  550 

_JJF_  (NELTL  , EO,  0)  GO  TO_  540  _ 

IF  (AtiS(ELTK)  ,LT,  1 .QE-mGO  TO  550 

GO  TO  Sap 

SPEC  IS  ZERO 

JP_L*BS( REMAIN)  .GT,  ,1E-1 1 )GO  TO  96 

IF  Sn'ELTL  .EG.  0)  GO  TO  565 

_IF(ABS(ELTK).GE.1.E-12)G0  TO  21 

IF(S1.LT.-1.0E-12)G0  TO  565 

GO  TO  545 

S2sELTK-Sl 

IF(S2)  5«5,5«5,?2  _ 

IF(ABS(S2/tLTK).LT,l,0E-l2)Gb  TO  545 

GO  TO  565  

SHUNTINTERVAL  IS  DUE  TO  ELT  BLOCK 

DELTs  SELT(nELT) 

TEX(«)si 

OXeOELT-  VAR(  1 ) 

REMAIN  s C I • D X 

REMAIKsREMAIN*Kl  

LTSHsl 

NELTsNELT»l 

NELTLsNELTL-1 

IF  (NELTL. t(J,0)G0  TO  560 

ELTKeKl*SELT(NELT) 

GO  TO  560  

DXsCI 

TEX  ( 3 ) s l 

GO  TO  560 

_IF  (NELTL  ,EU.  0)  GO  TO  98 
IF  (ELTK  ,LT,  (VAHKfRtMAIK) ) GO  TO  9« 
DXsREMAIN 
1EX(7)*1 

REMA I NSQ , Q 

GO  TO  560 

OX*CI _ 

TEx(3)*1 

_TEST*1  _ 

GO  TO  555 


BEGIN  RUNGE-KUTTA 


TEST*0 


I 


555  00  29  1 = 1, N 
_ 2«_&ULARClJ  =^yAH  (_Lu 

575  CUVAR (l)sVAR(l) 

576  00  25  1=1, N 

SOY ( I ) =DER  ( I ♦ 1 ) 


CUVAR ( 1 )=CUVAR( 1 )fDX/2,0 

CAtJL_  PE P SUB  .t 1 1 , CUVAR,  DE.R,  N,  V.AR1 

IF  (11  ,EU.  U)  GO  TO  120 

58Q_00_j>6_I  = 1 , N _ __ 

SOV(I)=SDY(I)*2,0*OER(I*1) 

»VAR  C I ♦ 1 )=S1  VA 


CALL  DERSUB  ( 1 1 , CUVAR, DER,N, VAR) 

IF  JL1 1 .EU  ,_«J_GU  JD _i iO_ 

585  DU  27  1=1, N 

_ _SPmj=SPY(JJ*2.0*UER{I*n 

27  CUV AR(Im  )=Sl  VAH(l)+Dx*OtR(  !♦!) 


UVUL 


BOUT] 


CALL  OERSUti  ( 1 1 , CU V A«, OER , N, V AR ) 

1E__L11_»E.Q  .__9J_Gj2_LQ_liii 

590  00  90  1=1, N 

___  J_PJL(11*15D  Ul  )90fcW(ltl)l/6.0 
90  CONTINUE 


nvi'iJWjnJK] 


28  DO30  1=1, N 

join  U SSO  Y U )_ 

YINCR(I)=0,0 

pep  uu . 

DXsOX/2,0 


IVI'Iil 


GO  TO  575 

TOUP  WAS  NOT  zerF 


29  DO  31  1*1, N 

NCR ( I ) =Y 1NCR ( 


IF  (LOOP  , EQ,  2)  GU  TO  33 

_ DO  32  1 = 1,  N _ _ 

SlVAR(I)=VAR(itl )+DX*YINCRlI) 

32  CUVAR(I*1 )=S1VAR(I) 

CUVAK(1)=VAR(1)+DX 
OOP* 


CALL  OERSUB  ( 1 1 , CUV AR , OEP, N, V AR ) 

IF_t.IJ_.tEG*  .9J_._GtL.I3  .120 

GO  TO  576 

33  LOUP=0  

H=2 • 0 *DX 

DO  39  1 = 1 ,N 

ERV0VH(I)  = ( YlNCR(I)/2,0-SDYl  (I)  )/15,0 

ERRVAL<I)*H*ERVOVM( Ij 

39  SlVAR(I)  = SlVAR(I)m*SDYCI  jURRVAL(I) 
C 

C S1VAR  HOLD  THE  approximate  ANSWERS 


IF  (SCIMAX)  36,35,36 
36  IF (ABS(SCIMAX-CI)  ,LTf 1,OE-12)GO  TO  38 
35  IF(ABS(H-CI),GT,1.0E-12)GO  TO  38 
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DCODE=0 

GO  TO  605  _ 

DCODfcM 

_C  ON  T 1NUE 

1*0 

!*!♦! 


IF  (I  , GT  , N)  GO  TO  a5 
IF_(ABS(S1VAR(I))  , t T , ELE2C I ) ) GO  T! 
RELER=AHS(FH»VAL(I)/S1VAW(I)) 

IF  (RELEF  ,GT.  SELF  HI))  GO  TO  615 
IF  (RELER  .GT.  RELMN(  1 ) ) DC0l>E  = 1 
GO  TO  ao 


CONTINUE 

IF  (OCUOE-1)  610,620,610 
CONTI NUt 

IF  (SSPEC)ai ,«2,ai  

IF  (SC  I V'A x ) a i , a j , a i 

Cls2.0+CI 

TEX (8)  = 1 

NHAF*NHAE-1 _____ 

GO  TO  620 

IF  (2,0*A8S(CI)  ,LEt_ABS(SCIMAX) ) GO 
CIsSCI^AX 

TEX ( 6 ) s 1 

GO  TO  620 


HALF  INTERVAL^ 

NHAFsNHAFfl  

TE  X (9 ) s 1 

NV ARs I ♦ 1 

IF  (NHAF-8)u7,a7,5Q5 
IF  (LTSM  ,EO,  0)  GO  TO  98 
TEST  = 1 

LTSH*0  _ 

NELT*NELT-1 

NElTLsNELIlM 


ELTKgKl*SELT (NtL  1 ) 
REMAINgO.O 

IF  (TPSH  , E 0 . 0)  GO  TO  U9 

TESTS  1 

TPsTP-SSPEC 

TRKsKl*TP 


TPSHsO 

IF  (SSPEC  , NE , 0.0)  GO  TO  996 
TESTsO 

IF  (ABS(C1-2.0*0X ) , GT,  t.E-1 
CIsDX 

DX=J)  X/2 , 0 

ClKcKHCl 

DO  50  I s 1 * N 

SIVARd  JsVAHCHl) 

DtR ( I ♦ 1 ) sF Ut R V ( I ) 

SOY  I (I)*YINC«(1)-SDV(I) 

Y1NCR(I)=0,0 

*K3s2 

IF  (ITEXT  .EQ.  1)  GO  TO  6J7 
LOOP* 1 


, GT,  1 ,E-12)GU  TO  1100 


TKX(12)sl 
501  KK3S4 

IF  ( 1 TEXT  ,EU.  lT  GO  TO  637 
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58  IF  (TEST  , t (J , 1)  GO  TO  520 

120  R£1U«  N 

C 

C RLdJMf  UlE  Interval 

C 


[ l'l'l 


OXsCI 
_.T  E a d 1 ) = 1 
k*3  = 3 


riimi 


nm«9M 


70  CIksCI*kj 

DO-fcO  1 = 1 __ 

DER ( I ♦ l JsFufcHVt 1 ) 

80 CUV AR  C I ) = VA«(I) 

CUVAfi(Ntl)s  VARir.+  l) 

IF  dPSn  . Ui *_J)l_GU_liULl 

TPsTP-SPEC 

J_RK8lP*Kl  

TPSHsO 

81_.IF  JLTSH  ,E0.  0)  GO  LQ  555  __ 
NELT*NElT-1 
Ms 


NELTL=NtLTL*l 

fclTKsSf  LI  (NELT)*K1  

GO  TO  555 

_6J6  >«  II  E ( 8 , 1 8 3 )_  V A H C 1 2.«_PX_  __ 

GO  TO  102 
IF ( Tfc  X ( 1 


1F(TEX(2).E0.U  WRITE  (6,172)  C 1 , C 1 MAX , SPEC 

6J7__JF<IEX(3).fcQ.i)  WRITE  (6,173)  

2F(TEX(«) .tQ. 1 ) WRITE  (6,17«)  M 

__ ILCI E x ( 5 ) . fc o • 1 ) wR  1 T t.  (6,175)  H 

IF (TtX(6) «E0. 1 ) WHITE  (6,176)  H 

1FCTEXC 7) .EO. 1 ) WRITE  (6. 184)  m 

IF(TEX(8).EO,l  ) WHITE  (6,177)  Cl 

IF(TEX(9).EU.I)  WHITE  (6,178)  N V A R , Cl ..  . 

IF(TEX(10).EQ.l)  wHITfc(6, 185)  NvAR,  OX 

UF_(  TEX(H).EO.l) WRITE(6,183)  VAR(l),  DX  

IF(TEX(12).EU.  1)  WRITE (6, 176)  vAR(  l) 

IF  ( TE XJJ37,E QJ] w RJ  T E ( b # 1 8 0 2 

Ic  f TEX (1 9 ) , tu,  1)  wH 1 Tt (6, 1 6 1 ) 

IF  t TEX ( 15) . EG. 1 ) wR  I ( E ( 6 , 182)  . 

102  DO  320  1=3,13 
320  TEX(I)so 

GO  TO  (120,99, 70,58) ,KK3 

171  FORMAT  (33"  INITIALIZATION  STARTS  A1  V AR  d2_=  f E 1 6 1 8/l_ 

172  FORMAT  (4H  C1*,E15.B,9H  C I M AX= , E 1 5 , 8 , 8H  SPEC=, E 1 5, 8/ ) 

1 73_F ORMAT  (37H  DX  ISJLmE  FULL  COMPUTING  INTERVAL^Cl/J. 

176  FORMAT  (28H  DX  IS  A SHORTENED  INTERVAL  ,tL5.6,25H  OuE  TO  A CRITIC 

1AL  VALUE/)  _ _ 

175  FORMAT  (28H  DX  IS  A SHORTENED  INTERVAL  #tl5,8,21H  DUE  TO  A SPEC  \ 
ALUE/ 


176  FORMAT  (28H  DX  IS  A S iORTENED  INTERVAL  ,E15.8,39H  DUE  TO  BOTH  A ; 

1PEC  and  critical  VALUE/)  

177  FORMAT  (27H  Cl  HAS  BEEN  LENGTHENED  TO  ,£16.6/) 
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TM  18 k 

178  FORMA T (5H  VAR(,I2,32H)  HAS  CAUSED  Cl  TU  BE  HALVED  TU  #Elb.6/) 

1 1?_FJ3RMAJ_ fiL7 h_VAR.( j.l_ ha S_M_E  N_UPD ATED  TO. El b ^Bx/J 

180  FORMAT  ( 3 1 H ERROR  RETURN-ELT  MOT  MGNOTONIC/) 

181_FQRMAT  (55H  ERROR  RETURN-HAyE  HALVED  9 TIMES  OVER  LAST  3 INTERV A L S 
1/) 

182  FORMAT  1 45H  ERROR  RETURN»HA vE  HALVED  9 CONSECUTIVE  TIMES/) 

183  FORMAT  ( 3 1 H INTERVAL  RECOMPUTED  AT  VAR  ( 1 ) = , E lt>.  8, 9h  rtltH  Oxs,El6,e 

1/1 

184  FORMAT (25H  DX  IS  SHORTENED  INTERVAL. E16.8.28H  DUE  TO  A PREVIOUS  EL 

__.„LL_mUt/L  .. 

18b  FORMAT (bH  VAR(,12,32H)  HAS  CAUSED  DX  TO  BE  HALVED  TO  .Elb.8,38H  BU 
IT  NOT  Cl  SINCE  Cl  ALREADY  SHORTENED/) 


500 


505 


510 


IERR=2 

TESTsO 
GOJLQ.  bl_ 
IERRc3 
TEXC15)sl 
TESTtO 

IERRS4 

= Q 

TEX(l«)sl 


SUBROUTINE  OEHSUB 

CGMhON/KUTTA/_  CUVAK(<?)  , VAR(2)  ,CI  ,0tR(2)  # II  

COMMON/F/RJtT , XNE, PE, AE,PINF  * AINF,UINF,RHINF 

COMMUN/0LAST/lBrt,Bi-.CaN#SREtO#«JtT2 

ROISTsCUVAR(2)/RJET2 

0ER(2)gAINF/SQRT  ( 1 . 0-SftCQN/ ( SQH TC ALOG ( ROI S T) ) *RD  I ST  ) ) 

RETURN 

END 


i 


) 
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